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I.  INTRODUCTION 


In  the  interest  of  building  an  axial  gain  cross  field  amplifier, 
analysis  was  desired  for  the  tapered  waveguide  (Fig.  1.).  A  stepped 
dielectric  profile  was  chosen  to  model  the  linearly  tapered  profile 
(Fig.  2).  The  formulation  presented  in  this  paper  is  based  on  the 
geometry  of  Fig.  3,  where  a  given  slab  is  removed  from  the  step.  Par¬ 
ticular  interest  in  the  behavior  of  certain  modes  in  the  presence  of  a 
discontinuous  dielectric  included  the  TEq  ^  TMq  ^  .  TEM,  and  several 

TEm,n  and  ™m,n  fflode8* 

The  original  approach  taken  for  the  problem  was  the  modal  expan¬ 
sion  technique.  This  traditional  approach  generates  a  dispersion  deter¬ 
minant  by  enforcing  continuity  of  certain  fields  across  an  interface 
between  two  mediums.  Roots  (kz  as  a  function  of  u>)  are  obtained  by 
finding  the  zeros  of  the  determinant  formed  from  these  boundary  condi¬ 
tions.  The  formulation  is  straightforward,  but  numerical  difficulties 
arise  when  it  is  translated  to  code  and  run.  This  motivated  the  need 
for  a  fresh  approach  to  the  problem. 

The  direct  integration  approach  simply  involves  integrating  two 
second  order  linear  differential  equations  for  TEm  n  and  T^  n  modes 
(m  *  0)  or  one  second  order  linear  differential  equation  for  TEq  n, 
TMq  n,  and  TEM  modes.  The  formulation  can  handle  discontinuous  and 
linearly  tapered  (as  a  function  of  r)  dielectric  profiles.  Theoreti- 


cally,  it 

can 

also  handle  combinations 

of  both 

to 

produce  a 

general 

dielectric 

profile. 

Irrelevant  of 

the 

number 

of 

dielectric 

layers 

comprising 

the 

profile, 

the  number 

and 

format 

of 

the  second  order 
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differential  equations  does  not  change  for  a  given  node.  There  are  no 
Bessel  functions,  and  consequently  no  need  to  evaluate  Bessel  function 
expansions  in  a  computer  program.  A  major  advantage  regards  the  ability 
to  analyze  the  TEM  mode  for  a  discontinuous  dielectric  profile,  as  shown 
in  Fig.  3A.  This  seems  logical  based  on  the  fact  that  all  modes  must 
obey  Maxwell's  equations  for  a  given  dielectric  profile.  In  stark 
contrast,  the  modal  expansion  technique  sheds  no  light  regarding  analy¬ 
sis  of  the  TEM  mode  for  a  discontinuous  dielectric  profile.  Finally, 
upon  completing  the  integration,  all  field  components  of  a  gi,ren  mode 
are  easily  computed  using  the  values  of  the  integration  parameters,  the 
mode  propagation  constant,  and  Maxwell's  equations. 

The  code  analyzes  the  geometry  of  Fig.  3A  and  its  effects  on  the 
following  mode  types:  TMjjj  n  (EH^  TEq  TMq  n,  and  TEM  modes.  The 
user  inputs  a  desired  dielectric  profile,  mode  type,  and  initial  condi¬ 
tions  from  which  the  propagation  constant  kz  and  all  relevant  field 
patterns  (as  a  function  of  r)  are  produced.  For  TEm  n  and  TMffl  n  modes, 
an  additional  parameter,  a,  is  produced  which  interprets  the  degree  of 
hybridization  due  to  the  dielectric  profile.  The  program  also  checks 
the  orthogonality  between  any  two  modes  of  a  given  mode  type  and  similar 
azimuthal  mode  number. 


TAPERED  COAXIAL  WAVEGUIDE 


INNER  CONDUCTOR 

DIELECTRIC 

AIR 


OUTER  CONDUCTOR 


Fig.  1.  Tapered  coaxial  waveguide 
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STEP  APPROXIMATION 
OF  THE 

TAPERED  COAXIAL  WAVEGUIDE 


TYPICAL  SECTION  USED  IN  FORMULATION 


Fig*  2.  Step  approximation  of  the  tapered  coaxial  vavegulde. 
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MODAL  EXPANSION  APPROACH 


The  modal  expansion  technique  Is  presented  In  this  chapter. 
Section  2.1  derives  the  dispersion  determinant  for  the  geometry  of  Fig. 
3A.  We  start  with  the  scalar  Helmholtz  equation  and  proceed  to  separate 
the  variables.  The  scalar  wave  functions  (eigen  functions)  are  con¬ 
structed,  which  then  allows  the  appropriate  boundary  conditions  to  be 
applied  at  the  dielectric  Interface  (r  ■  b).  Finally,  the  dispersion 
determinant  is  found  which  describes  the  modes  propagating  in  our  sys¬ 
tem.  Section  2.2  discusses  the  drawbacks  using  this  approach.  We  note 
that  primes  in  the  equations  represent  the  derivative  with  respect  to 
the  radial  coordinate  r. 

2.1  Derivation  of  the  Dispersion  Determinant 

The  derivation  that  follows  is  based  on  that  of  Harrington1  and 
Rothwell.2  The  scalar  wave  functions  for  the  geometry  of  Fig.  3  must 
obey  the  scalar  Helmholtz  equation  written  here  in  cylindrical  coordi¬ 
nates: 

1  3_ 
r  3r 

where 

'll  ■  scalar  wave  function 


jhp 

3r 


1  32i|» 

7  77 


32<1> 


k2i|» 


(1) 


3z 


k  ■  constant 

The  common  approach  taken  to  solve  Eq.  1  is  to  separate  the  variables. 
The  usual  form  for  the  solution  is 


-  6 


-  R(r)#($)Z(z) 


Following  the  standard  procedure,  the  separated  equations  are 


h  [r  g]  *  [V  *  "2] 


i-t  -  „2.  -  0 

d<T 


—  +  k2Z  *  0 
.  2  z 
dz 


where 


2  2  2 

kz  +  lc  -  k 
z  r 


Solutions  to  Eqs.  4  and  5  are  harmonic  functions.  Equation  3  is 
known  as  Bessel's  equation  of  order  m  and  argument  k^  •  r  with  solutions 
having  the  general  form 


R(r)  *  Z  (k  ,  r.  A,  B)  2  AJ  (k  r)  +  BY  (k  r) 
or  r’  J  or  r  '  nr  r  J 


where 


J  (k  rl  *  Bessel's  function  of  the  first  kind 
nr  r  1 


Y  fk  r)  *  Bessel's  function  of  the  second  kind 
nr  r  1 


A,  B  *  constants 
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and  kr  may  be  real  or  imaginary, 
form 

*(♦)  - 
*(*)  " 
Z(z)  - 


Solutions  to  Eqs.  4  and  5  are  of  the 


cos  (m$) 

(7) 

sin  (m$) 

(8) 

-jk  z 

z 

e 

(9) 

where  k  is  the  propagation  constant.  The  scalar  wave  functions  for 
z 

each  region  of  Fig.  3A  are: 


Region  1 

” jk  z 

-  Z^nl^(k  .  r,  A  ,  A.)  cos  (m<t>)  e‘  Z  (10a) 

m  m  ri  l  z 

.  i .  .  .  -jk  a 

*m  -  ZoJ(kri*  B2)  sin  (m$)  e  z  (10b) 

Region  2 

-  Zlm2^kr2’  r’  Cl  »  C2 )  cos  e  2  (Ha) 

*  jk  z 

<f^e2^  -  Zme2^kr2’  r>  Dl*  °2^  sin  6  *  (Hb) 


where  A^»  B^,  and  are  constants  (i  -  1,  2).  The  superscripts  ml  and 
m2  denote  the  contributions  from  TM  modes  for  regions  1  and  2;  the 


superscripts  el  and  e2  denote  contributions  from  TE  modes  for  regions  1 
and  2.  The  separation  equations  which  the  wave  functions  must  satisfy 


(12a) 


z 


“2Vi 


z 


k2  "  “2V2 


(12b) 


for  regions  1  and  2,  respectively.  The  resultant  electric  and  magnetic 
fields  in  each  region  have  contributions  from  both  the  TE  and  TM 
modes.  Consequently,  hybrid  modes  will  propagate  in  our  system. 

Boundary  conditions  were  applied  to  the  following  field  components 
written  in  general  for  region  i(i  *  1,  2): 


%  - 


5?  2lmI)(krl-  r’  "l-  *2)  +  kriZ;(el>(kri-  r-  *3'  V 


sin  (m$)  e 


-jk.z 


(13a) 


k2  , 

m  _£i_  z(mi)(kri’  r*  K5’  K6'  cos  e 
z  juEj  m  ^ 


(13b) 
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k  Z,(ml) 
rl  m 


(k 


ri* 


*7 


mk 

,  K  )  +  -5-  Z 

8  wu^r 


(ei) 

m 


(krl.  V  t10) 


-jkzz 

cos  (tn$)  e 


(13c) 


H 

z 


^-Z(el)(k  . 

j  tou1  m  ri 


-jk  z 

r,  Kn.  K12)  sin  (m«t> )  e  2 


(13d) 


where  1^  (i  ■  1,  ...,  12)  is  a  constant. 
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constants 


at  r  •  a  and 


(17a) 


n  ! 2 
S' 

1  b8 


(17b) 


at  r  *  c.  Using  Eq.  6,  we  eliminate  A2,  B2,  C2,  and  D2  with  the  aid  of 
Eqs.  16  and  17, 


}  "  AlZm^krl’  r’  lf  _S1S2  ^ 


(18a) 


♦I"0  "  BlZm^krl*  r*  ‘Si(S2^'^ 


(18b) 


for  region  1.  and 


ClZm(kr2’  r’  _S7S8  ^ 


(19a) 


*i°2)  ■  DlZ»Ckr2-  r-  'S7(S8^‘J 


(19b) 


for  region  2.  At  the  dielectric  interface  (r  ■  b)  E^,  Ez,  H^t  and  Hz 
must  be  continuous.  Using  equation  sets  of  Eqs.  13  and  14  for  each 
region  produces: 


-  11  - 


m 


XV’V-VJVUWrfVjf 


t  <  I  .  /  *’-**• 


E  Continuous 
z 


H  Continuous 


To  further  simplify  Eqs.  20  through  23,  we  define  the  following 


stants 


With  the  definition  of  Eq.  24,  Eqs.  20  through  23  become 


mCjkjj 


G4krl 


mG5kz 


G8kr2 


A  _ _  +  B  -  C  --  -  -  D  — — — 

1  <*>£^82  1  S2  1  i^bSg  1  Sg 


G  k2  G  k2 

A.  -  C  -Lli  -  0 

1  eis2  1  e2s8 


G3krl 


nC2kz 


G6kr2 


nG7kz 


A  — - —  -  B  - —  -  +  C  — -r —  +  D  ■  -  -  -  ■  0 

1  S2  1  ujyjbS^  1  Sg  1  u^bSg 


G.k2,  G  k2 

b  -2y-1-  -  D  “  0 
1  U1S2  1  U2S8 


for  E^,  Ez ,  H^,  and  Hz,  respectively.  Based  on  Eqs.  25  through  28,  the 
characteristic  equation  in  determlnantal  form  is 
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.  f*  v 
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°2krl 


°7k?2 
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G5kr2 


G3krl 


mG„k 
2  z 


G6kr2 


mG-,k 
7  z 


toe  bb 
1  2 


toe2bS8 


This  determinant  represents  solutions  for  hybrid  modes  and  columns 


marked  TM  or  TE  denote  the  relative  contributions  to  the  hybridiza¬ 


tion.  In  general,  the  zeros  of  this  determinant  lie  in  the  complex 


plane. 


If  we  set  m  -  0,  the  determinant  simplifies  to 


G2krl 


°7kr2 


Glkrl 


G5kr2 


-  0  (30) 


c3kri 


G6kr2 


G4krl 


G8kr2 


Determinant  Eq.  30  is  the  uncoupled  version  of  Eq.  29  and  it  represents 
solutions  to  pure  TE  and  TM  modes.  This  is  illustrated  by  interchanging 
the  first  and  third  rows  of  Eq.  30,  producing 


where 


for  TM  modes,  and 


for  TE  modes. 


°3krl 


Glkrl 
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°2krl 


G.k  , 
4  rl 


G^k  o 
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G„k 

8  r2 
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2.2  Drawbacks 


As  dielectric  layers  are  added  to  the  profile,  the  increasing 
number  of  continuity  conditions  applied  at  each  interface  results  in  a 
dispersion  determinant  that  grows  as  the  square  of  the  number  of  transi¬ 
tions.  This  drawback  has  two  consequences.  The  first  is  that  the 
dispersion  determinant  must  be  rederived  for  each  new  layer  of  dielec¬ 
tric  added.  The  second  is,  obviously,  the  increasingly  complex  deter¬ 
minant  for  which  a  zero  must  be  found. 

The  modal  expansion  approach  is  a  technique  where  continuity  of 
certain  field  components  is  applied  at  the  boundary  where  the  character¬ 
istics  of  the  medium  change  abruptly.  Cases  where  the  dielectric  is 
tapered  continuously  as  a  function  of  r  cannot  be  handled  by  this  formu¬ 
lation.  Another  drawback  results  from  numerically  approximating 
Bessel's  function  in  a  computer  program.  Problems  arise  when  the  argu¬ 
ments  of  the  expansions  used  approach  zero  due  to  kr  approaching  zero. 
Instabilities  and/or  numerical  overflow  result,  placing  restrictions  on 
investigating  modes  which  are  in  transition  between  the  fast  and  slow 
wave  regions  of  the  structure.  Due  to  the  above  mentioned  problems,  we 
could  not  analyze  the  dielectric  loaded  coaxial  waveguide  success¬ 
fully.  This  led  us  to  develop  the  technique  described  in  the  next 


chapter. 


III.  DIRECT  INTEGRATION  OF  MAXWELL'S  EQUATIONS 


This  chapter  will  lay  the  foundations  used  to  calculate  the  elec¬ 
tromagnetic  dispersion  of  a  coaxial  waveguide  with  an  arbitrary  radial 
dielectric  profile.  Starting  with  Maxwell’s  equations,3  Section  3.1 
will  develop  a  coupled  second  order  differential  equation  system  in  the 
radial  coordinate.  The  propagation  constant  sought,  kz,  is  found  by 
integrating  this  system  of  differential  equations  at  a  given  frequency 
(to)  and  azimuthal  eigenvalue  (m).  Integration  is  performed  using  a 
"shooting"  method,  where  the  shooting  parameter,  kz,  is  varied  until 
certain  boundary  conditions  are  satisfield.  Section  3.2  discusses  the 
initial  and  boundary  conditions  used  in  the  integration.  For  discontin¬ 
uous  dielectric  profiles.  Section  3.3  formulates  the  jump  conditions  for 
the  equations.  The  equations  for  the  field  components,  which  are 
directly  obtained  from  the  solutions  to  the  integration,  are  given  in 
Section  3.4.  The  normalization  of  the  field  components  and  computation 
of  orthogonality  are  discussed  in  Section  3.5.  Finally,  in  Section  3.6, 
the  equations  are  transformed  into  a  format  suitable  for  numerical 
integration.  We  note  that  boldface  type  will  represent  vector  quanti¬ 
ties. 


3 . 1  Derivation  of  the  Differential  Equation  Systems 

The  general  Maxwell's  equations  for  a  region  free  of  charges  and 
currents  (p  *  0,  J  =  0)  are 


7  •  E  «  0 


(31a) 
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7  •  B  -  0 


(32a) 


7  x  H 


3D 

at 


7  X  B 


3B 

at 


For  linear  dielectric  media,  the  following  relation  holds: 


(33a) 

(34a) 


D  -  cK  (35) 

The  formulation  will  be  limited  to  the  case  where  u  »  Uq  (w^  is  the 
permeability  of  vacuum): 


B  -  uqH  (36) 

Assuming  harmonic  variation  and  uniform  propagation  in  the  +z  direc- 
j(wt-kzz) 

tion,  e  ,  Gqs.  31a  through  34a  become 


7  •  D  -  0  (31b) 
7  •  H  -  0  (32b) 
7  x  H  -  juD  (33b) 
7  x  K  -  -juw0H  (34b) 


where  Eq.  36  was  used  in  Eq.  34b  for  B.  Rewriting  Eqs.  31b  through  34b 
in  cylindrical  coordinates  yields 


3D  ,  3D. 

_ E  +  I _ 1  + 

3r  r  3<> 


0 


(31c) 
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H  (r ,  <fr)  ♦  H  (r)  cos  (tu<f>) 
9  9 


(37h) 


3 


H  (r,  $)  -*•  H  (r)  sin  (m$) 
z  z 


(371) 


Using  this  polarization  in  Eqs.  31c  through  34c-z  and  canceling  common 
cosine  and  sine  terms  gives 


3D_  mD.  D_ 


(31d) 


3H_  raH.  H_ 


(3  2d) 


-  +  jk  H  -  juD  -  0 

r  J  z  $  J  r 


-7 - jk  H  -  ju»D.  -  0 

3r  J  z  r  J 


(33d-r) 


(33d-<fr) 


3H.  mH..  H. 

**.-o 


(33d-z) 


—  +  jkzE+  +  ja,u0Hr  -  0 


(34d-r) 


JT  -  jkzEr  +  ^0H*  "  ° 


Ei 

+  — -  +  -i  +  jup.H  -  0 
3r  r  r  J  0  z 


(  34d— 4> ) 


( 34d-z) 
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mm 


Equations  31d  through  34d-z  will  provide  the  basis  for  determining  a 
system  of  second  order  linear  differential  equations.  Before  proceed¬ 
ing,  we  further  assume  that  e  Is  a  function  of  the  r  coordinate  only. 

3Ez 

We  begin  by  solving  Eq.  34c-$  for  and  then  using  the  constltu- 

3Dz 

tive  relation  (Eq.  35)  to  generate  which  yields 


3r 


-  Ez  If  -  J-VS  *  JkzDz 


(38) 


If  we  solve  Eq.  31d  for  jDz  and  then  take  the  radial  derivative,  we  have 


3(J°2) 

5r 


32D  ,  3D  D  3E 

r  ,  1  r  r  .  me  d> 

7 T"  +  7  Tr  T  +  ~  Tr“ 
3r  r 


mE ,  „  meE, 
+  _ _ 4» 


W 


1 

F" 
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(39) 


where  the  constitutive  relationship  for  D.  was  used.  By  taking  Eq.  39 

3Dz  2  2 

and  substituting  for  in  Eq.  38,  we  can  solve  for  3  D  /3r  , 


32D_ 


3r 


,  3D  D  3E.  mE.  meE.  _ 

1  r  r  me  9  6  3e  ®  _  3e  ,  „  ,  2_ 

—  - —  +  — - - -  -  — -  ^  +  jk  E  - my  ek  H  +  k  D 

r  3r  r2  r  3r  r  3r  f2  z  z  3r  0  z  $  z  r 

(40) 


In  order  to  put  Eq.  40  Into  a  form  involving  only  electric  field  com¬ 
ponents  and  their  derivatives,  we  substitute  for  and  Hz  from  Eqs. 
33d-r  and  34d-z,  respectively, 


32D 


3r 


,2  .  ,2  m2  +  11  3 ( lne ) 

3D 

+  _ - 

3  (lne)  1 

2meE . 
+  * 

kz  k  +  2  r  3r 

3r 

3r  r 

+  j- 

r 

r 

(41) 


where  the  following  relations  were  used: 
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Equation  41  represents  the  first  half  of  our  differential  equation 


system. 

To  generate  the  other  half  of  our  differential  equation  system,  we 

2  2 

take  the  radial  derivative  of  Eq.  34d-z  and  solve  for  3  E  /3r  : 
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3D  mD_ 


w  m  r  ^  r  1  3e  1  “4  3e 

72  "  “7l3F“  +  1""T37Tr-7Tr"7t3F+2 

3r  r  e  r 


-  jmy 


3Hz  ^0Hz  3, 


0  3r 


where  the  constitutive  relation  for  Er  was  used.  Using  Eqs.  33d-<fr  and 

3Hz 

34d-z,  a  solution  for  -sr—  as  a.  function  of  E,  and  Ez  is  obtained. 
3Hz 

Substituting  -jj—  into  Eq.  44  yields 
2  _ 

!_%  m  E  k2  _21  _  1  3(lne)  _  ^(lne)  l] 

~ T  "  E*  kz  k  +  2  r  3r  3r  3r  +  rj 

3r  L  r  J  L  J 


_  3D_  mD_  jmk D 

m  r  .  r  .  J  z  a 


re  3r 
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where  Eqs.  42  and  43  were  again  used.  Equation  45  is  finally  put  into  a 
form  involving  only  electric  field  components  by  substituting  for  Dz  and 
H  from  Eqs.  31d  and  34d-z,  respectively: 


7T  "  e4  kz  "  + 


m2  +  1  _  1  3E»  ^r  [  2  3(  lne ) 

2  r  3r  2  r  3r 

r  J  r  e  L  J 


Equation  46  completes  the  formulation  of  the  differential  equation 
system. 


Upon  normalizing  Eq.  41  by  dividing  by  (permittivity  of  vacuum) 
in  order  to  have  consistent  units  with  Eq.  46,  the  second  order  differ¬ 
ential  equation  system  is 


3r2  C0 


Dr  ,  2  ,  2  .  m2  +  1  .  1  3(lne) 

ert  Kz  ~  *  2  r  3r 

0  r 
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(47a) 
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2  3E,  mD  P_  ...  . 

m  +  1  J_  4  _ 2^  3(lne) 

2  "  r  3r  re  *r  3r 

r  J  L  J 


(47b) 


where  e  -  (e^  is  the  relative  dielectric  permittivity).  Equation 
47  is  used  to  find  the  propagation  constant,  kz,  for  the  TEm>n  and  TM^  n 
modes  as  well  as  the  "TEm>n  like"  (HEm>n)  and  like"  (EH^  n)modes. 
For  the  special  case  of  m  ■  0,  Eq.  47  becomes 


2  .  2 


1  r  3(lne)  1 


Equation  48  is  the  uncoupled  form  of  Eq.  47,  with  Eq.  48a  applying  to 
TMQ^n  and  TEM  modes,  and  Eq.  48b  applying  to  TEgfn  modes.  It  is  inter¬ 
esting  to  note  the  analogous  relationship  between  Eqs.  47  and  48  and  the 


dispersion  determinants  of  Eq.  29  and  30. 


3. 2  Initial  and  Boundary  Conditions 


For  TEg^n,  TMg^n,  anc*  TEM  modes,  a  one-dimensional  "shooting" 
method  is  used,  with  k2  being  the  shooting  parameter.  However,  "TEmjtl 
like"  and  “TMm^n  like"  modes  require  a  two-dimensional  "shooting” 
method.  The  shooting  parameters  here  are  kz  and  the  dimensionless 
parameter  a,  which  we  define  as 


a  = 


JH, 


C0°r 


(49) 


at  r  -  a  and  eg  *  speed  of  light  in  vacuum.  The  degree  to  which  a  mode 
has  characteristic  TE  behavior  (H^  component  dominant)  or  characteristic 
TM  behavior  (E^  component  dominant)  is  represented  by  a.  The  initial 
magnitude  of  a  points  the  integration  in  the  direction  of  an  HEm>n  or 
EHm,n  mode  when  a  nonuniform  dielectric  is  present.  When  the  boundary 
conditions  are  satisfied,  the  final  magnitude  of  a  describes  the  degree 
of  actual  hybridization.  There  are  two  limiting  cases  for  a  which  occur 
when  a  nonuniform  dielectric  profile  approaches  a  uniform  profile: 


I 

i 


a|  ♦  0 


as  pure  TM  modes  are  approached,  and 


a|  ~  10 


as  pure  TE  modes  are-  approached. 
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Integration  of  Eqs.  47  and  48  requires  the  value  of  a  given  field 
component  and  its  derivative  with  respect  to  the  radial  coordinate  at 
the  inner  conductor  of  Fig.  3A  (r  ■  a).  The  tangential  electric  field 


components  must  satisfy  the  boundary  condition 


E  -  0  (50a) 

$ 


E  -  0 
z 


(50b) 


at  r  *  a  and  r  »  c,  respectively.  Using  this  result  in  Eq.  31d  and  34d- 
z  yields 
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(51) 


3E* 

3r 
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(52) 


Using  Eqs.  49  and  the  constitutive  relation  for  Er,  Eq.  52  becomes 


(53) 


Summarizing  the  initial  conditions  for  uniform  and  nonuniform  dielectric 
profiles: 

1.  Equations  50b  and  51  are  used  for  the  TMg^n  and  TEM  modes. 

2.  Equation  50a  is  used  for  TEq  n  rodes. 

3.  Equations  50,  51,  and  53  are  used  for  TEm  n  (HEn  n)  and 
™m,n  (EH^J  modes. 
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As  mentioned,  the  "shooting'*  method  varies  the  shooting  parameters 
until  the  boundary  condition  at  r  *  c  is  satisfied.  For  TM()jn  an<* 
modes,  Eq.  50b  must  be  satisfied,  while  TEq  n  modes  require  Eq.  50a  to 
be  satisfied.  Finally,  Eqs.  50a  and  50b  are  combined  to  form  the  bound¬ 
ary  condition  for  TEmjn  (HEm>n)  and  TM^  (EH  )  modes: 

-  0  (54) 


We  note  that  for  the  two-dimensional  shooting,  kz  enters  only  in  the 
differential  equations  and  does  not  appear  in  the  initial  or  boundary 
conditions.  Conversely,  a  enters  only  in  the  initial  conditions  and 
does  not  appear  in  the  differential  equations. 


3.3  Discontinuous  Dielectric  Profile  Formulation 


where  superscripts  1  and  2  denote  the  dielectric  regions.  Substituting 

for  and  E^,  using  Eq.  31d,  we  generate  Eq.  56  from  Eq.  55c, 
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where  e^  and  e^  represent  the  dielectric  permittivities  of  regions  1  and 

[21 

2,  respectively.  Solving  Eq.  56  for  3D^  /3r, 
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Similarly,  we  substitute  for  and  H^,  using  Eq.  34d-z  into  Eq. 


Solving  Eq.  58  for  — ji — f 

„[21 


\C2  ' 
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For  the  special  case  of  u  ■  0,  Eq.  59  becomes 
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Summarizing  Eqs.  57,  59,  and  60  for  reference  (Eq.  57  is  divided  by  e^). 
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Equation  61  is  used  for  TEra>n  (HEffl  n)  and  TMm>n  (EH^  n)  modes. 
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Equation  62  is  used  for  TMgfn  and  TEM  modes. 
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Equation  63  is  used  for  TEg>n  modes. 

3.4  Field  Equations 

The  following  field  components  may  be  computed  once  the  param- 
Dr  1  3Dr  3E* 

eters  — ,  E. »  —  -gj—  ,  and  have  been  determined.  Equations  31d, 
0  9  0 

34d-z,  34d-r,  and  33d-r  are  used  to  derive  Eqs.  64  through  67,  respec¬ 
tively, 
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3.5  Normalization  and  Orthogonality  Relations 


Once  computed,  all  field  components  are  normalized  with  respect  to 
Poynting’s  vector  taken  over  the  cross-sectional  area  of  the  waveguide: 
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where 


6.  ,  *  Kronecker  delta  function 
i,k 

dA  -  r  dr  d$ 

A 

n  *  unit  vector  normal  to  dA  in  +z  direction 


-  29  - 


Normalization  for  a  given  mode  requires  computing  Eq.  68  with  £  m  k 


using  the  unnormalized  fields,  taking  the  square  root  of  the  absolute 
value  and  dividing  each  component  by  the  result. 

Orthogonality  between  two  given  modes  £  and  k  (£  *  k)  is  an  impor¬ 
tant  condition  that  must  be  satisfied  in  order  to  have  any  confidence  in 
the  solutions  for  kz  ^  and  kz^.  Equation  68  is  used  to  evaluate  the 
orthogonality  between  the  normalized  modes  £  and  k.  Due  to  the  normali¬ 
zation,  Eq.  68  gives  the  value  of  1  when  self  orthogonality  is  computer 


(£  -  k). 


3.6  Form  for  Numerical  Integration 


Integration  of  Eqs.  47  and  48  is  facilitated  if  they  are  trans¬ 
formed  into  a  dimensionless  form.  We  define  the  following  normaliza¬ 
tion: 


r  «  — 
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where 


a  <  r  <  c 


As  a  consequence  of  Eq.  69a,  we  have 
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Using  Eq.  69,  Eq.  47  becomes 
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Similarly  for  Eq.  48,  we  have 
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Equations  70  and  71  are  dimensionless. 

The  jump  discontinuity  equations,  Eq.  61,  62,  and  63,  are  also  put 
into  dimensionless  form,  since  they  are  used  in  the  integration  process: 
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The  IMSL  routine  DGEAR  (see  Appendix  A.l)  is  used  to  integrate 


Eqs.  70  and  71.  But  DGEAR  requires  a  first  order  linear  differential 


equation  system.  To  accomodate  DGEAR,  Eqs.  70  and  71  must  be  trans¬ 


formed  into  a  first  order  linear  differential  equation  system. 


For  TEn,>n  (HEm>n)  and  TM^n  (EHffl  n)  modes,  we  define  the  following 


variables: 
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where  the  prime  denotes  the  derivative  with  respect  to  the  normalized 


radial  coordinate.  With  the  aid  of  Eqs.  75  and  76,  Eq.  70  becomes 
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(77a) 


(77b) 


(77c) 


(77d) 


Equation  77  constitutes  the  first  order  differential  equation  system 
which  is  integrated  by  DGEAR. 

Similarly,  for  TMg  n  and  TEM  modes,  we  define 


With  the  definition  of  Eqs.  78  and  79,  Eq.  71a  is  transformed 
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Finally,  for  TEoftl  modes,  we  define 
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Transforming  Eq.  71b  using  Eqs.  81  and  82, 
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Equations  80  and  83  represent  first  order  linear  differential  equation 
systems  integrated  by  DGEAR. 

Another  requirement  for  DGEAR  is  a  system  of  partial  derivatives 

P®i , J  defined  as  the  partial  derivative  of  YJ  with  respect  to  Yj.  In 
light  of  this  definition,  Eq.  77  generates 
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Equation  84  represents  the  partial  derivative  system  for  TEm  n  (HEm  n 
and  TMrafn  (HEm  n)  modes.  Repeating  the  process  for  TMqj0  and  TEM  modes 


Eq.  80  produces  the  partial  derivative  system. 
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Similarly  for  TEg  modes,  Eq.  83  produces  the  partial  derivative  sys- 
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With  the  first  order  differential  equation  and  partial  derivative  equa¬ 
tion  systems  constructed,  we  now  turn  our  attention  to  initial  and 
boundary  conditions.  In  each  case  to  follow,  we  are  free  to  choose  the 
value  of  one  variable,  since  it  affects  only  the  magnitude  of  the 
results* 

Dr 

Choosing  —  ■  1  coupled  with  the  aid  of  Eqs.  50,  51,  and  53,  Eq. 
0 

75  yields  the  initial  conditions  for  TE„  _  (he  1  and  TM_  _  (EH  ) 
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where  e  is  the  relative  permittivity  at  the  boundary  (Fig.  3).  The 
r 

corresponding  unnormalized  field  components  are 
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For  TMq  and  TEM  modes,  we  again  choose  —  ■  1  and  use  Eqs.  50b  and  51 
’  0 
from  which  Eq.  78  yields 
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The  corresponding  unnormalized  field  components  are 


(90a) 


Finally,  for  TEg>n  modes, 
which  Eq.  81  yields 
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1  and  the  use  of  Eq.  50a  from 
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The  corresponding  unnormalized  field  components  are 
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The  boundary  conditions  that  must  be  satisfied  at  r  *  c  in  the 
shooting  method  are  taken  from  Eqs.  50  and  51.  For  TMgfn  an<*  TEM  modes, 
we  have  from  Eq.  51, 
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Equation  50a  is  used  for  TEo,n  ra°des, 
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Finally ,  for  TEm>n  (HEm>n)  and  TMm,n  (EH^J  modes,  Eq.  54  yields  the 


boundary  condition 
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IV.  PROGRAM  OUTLINE  AND  USE 


Chapter  4  presents  the  key  aspects  of  the  computer  program  used  to 
implement  the  theory  and  equations  of  Chapter  3.  Section  4.1  starts 
with  an  overview  of  the  program,  showing  the  general  flow  of  logic.  The 
techniques  used  to  handle  a  general  dielectric  profile  are  discussed  in 
Section  4.2.  The  process  for  finding  a  mode  is  presented  in  Section 
4.3. 


4.1  Program  Flow 

The  Fortran  code  is  laid  out  in  blocks  with  a  number  and  comment 
delimiting  each  block.  If  necessary,  a  block  will  include  a  brief 
explanation  of  its  function.  The  blocks  making  up  each  procedure  (main 
program  and  subroutines)  are  presented  below.  Any  variable  used  in  the 
discussion  is  defined  in  Appendix  D. 

Main  Program  (C0AD7R) 

1.00  INITIALIZATION  SECTION 

Initialize  constants  and  logicals  used  throughout  program. 
1.10  VARIABLES  USED  IN  ZREAL1 

Initialize  variables  used  in  the  calling  argument  of  the 
zero  finding  routine  ZREAL1. 

1.15  VARIABLES  USED  IN  E04JBF 

Initialize  variables  used  in  the  calling  argument  of  the 
minimization  routine  E04JBF. 
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VARIABLES  USED  IN  DCADRE 


Initialize  variables  used  in  the  calling  argument  of  the 
integration  routine  DCADRE. 

SET  DEFAULTS 

Initialize  constants  and  logicals  which  are  input  parame¬ 
ters. 

OPEN  NAMELIST  AND  PRINT  FILES 

Input  variables  are  read  in  via  the  NAMELIST  file.  The 
print  file  C0AD7P  is  opened. 

CHECK  MODE  LOGIC 

A  check  is  made  to  ensure  that  only  one  mode  type  is 
chosen. 

DIELECTRIC  TAPER 

If  desired,  a  linearly  tapered  dielectric  profile  is  con¬ 
structed. 

INITIALIZE  DIELECTRIC  AND  RADIAL  ARRAYS 
DEFAULT  FOR  NJUMP  EQUAL  TO  ZERO 

Initialize  NJUMP  equal  to  zero.  Set  default  values  for 
DRATIO(l),  AJUMP(l)  and  IJUMP(l). 

EXAMINE  DIELECTRIC  PROFILE  FOR  ANY  DISCONTINUITIES 
Assign  appropriate  values  to  pertinent  arrays  when  a  dis¬ 
continuity  is  encountered. 

PRINT  OUT  PHYSICAL  CONDITIONS 


1.60  SMOOTH  DRELD  AND  DRELN 


Discontinuities,  if  any,  in  DRELD  and  DRELN  are  removed. 
1.65  INITIALIZE  DIFFERENCE  ARRAYS 

1.70  GENERATE  SPLINE  COEFFICIENTS  FOR  DRELA  AND  DRELN 
1.90  OUTER  LOOP  ON  MODE  NUMBER 

This  loop  evaluates  K  distinct  modes. 

2.00  OUTER  FREQUENCY  LOOP 

This  loop  finds  the  propagation  constant  kz  at  L  distinct 
frequencies  for  the  kth  mode. 

2.10  CASE  FOR  TEQ>n,  TMq  ,  OR  TEM  MODES 

The  propagation  constant  k  for  the  TEn  ,  TMn  ,  or  TEM 

z  u  j  n  u  j  n 

mode  is  sought. 

2.15  CALL  ZREAL1 

The  zero  finding  routine  ZREAL1  is  called. 

2.20  CASE  FOR  TEm  n  OR  TMm  mode  is  sought. 

2.30  CALL  E04JBF 

Minimization  routine  E04JBF  is  called. 

4.00  EVALUATE  FIELD  COMPONENTS 

Depending  on  the  mode  type  chosen,  initial  conditions  are 
set  for  one  of  the  following  pairs  of  equations: 

TEjj,  n  and  TM^  n  modes  :  equations  87  and  88 

TMQ^n  and  TEM  modes  :  equations  89  and  90 

TEq  n  modes  :  equations  91  and  92. 

4.10  SET  UP  FOR  DGEAR 

Initialize  variables  used  in  the  calling  argument  of  the 
integration  routine  DGEAR. 


SET  UP  FOR  NEXT  INTEGRATION 

Depending  on  the  mode  type  chosen  and  if  integration  is 
currently  at  a  discontinuity  in  the  dielectric  profile, 
compute  one  of  the  following  cases: 


TEm,n 

and 

TM_  „  modes 
m,n 

:  equations 

61 

and 

88 

™0,n 

and 

TEM  modes 

:  equations 

62 

and 

80 

TE0,n 

modes 

:  equations 

63 

and 

82 

UNNORMALIZED  FIELDS 

Depending  on  the  mode  type  chosen,  compute  the  remaining 
field  components. 

Define  difference  field  arrays  for  one  of  the  following 
mode  types: 

TEm  n  and  TMm  n  modes  :  ERD(l),  DTD(l)  ,  DZD(l) 

TM0n  and  TEM  modes  :  ERD(l),  DZD(l) 

TEq  n  modes  :  DTD(l) 

SMOOTH  DIFFERENCE  FIELDS  FOR  SPLINE  PURPOSES 
If  necessary,  the  difference  field  arrays  will  be 
"smoothed"  for  one  of  the  following  modes: 

TEm  n  and  ™m  n  modes  :  ERD(l),  DTD(l),  DZD(l) 

™0,n  and  TEM  mode8  ;  ERD(l),  DZD(l) 

TEq  n  modes  :  DTD(l). 

CALL  NORMALIZING  ROUTINE 


Normalize  field  components  of  a  chosen  mode  type. 


7.00  END  OF  OUTER  FREQUENCY  LOOP* 

7.20  PRINTOUT  RESULTS 

For  a  given  mode  type,  write  appropriate  results  to  file 
COAD7P. 

7.40  PLOT  DATA 

For  a  given  mode  type,  plot  appropriate  graphs  of  field 
components  and  write  to  file  C0AD7P. 

7.60  CALL  GRAPHING  ROUTINE 

For  a  given  mode  type,  appropriate  field  components  are 
plotted  and  hard  copy  printouts  are  made. 

7.90  END  OF  OUTER  LOOP  ON  MODE  NUMBER 
7.95  ORTHOGONALITY  CHECK 

Orthogonality  between  two  chosen  modes  is  evaluated  by 
computing  Poynting's  vector  over  the  cross-sectional  area. 
8.00  TERMINATION 
9.00  FORMAT  STATEMENTS 

Subroutine  DERIV 

The  integration  routine  DGEAR  calls  subroutine  DERIV,  which 
defines  and  evaluates  the  derivatives  of  the  second  order  linear  differ¬ 
ential  equation  system  (Eqs.  47  and  48). 

2.00  DEFINE  DIELECTRIC  AND  FIRST  DERIVATIVE  VALUES 

This  section  evaluates  FDRV  and  DREL  based  on  the  present 
radial  position  used  by  the  integration  routine  DGEAR. 


2.20  EVALUATE  BERIVATIVES 

Compute  one  of  the  following  first  order  differential 
equation  systems  for  a  given  mode  type: 

TE^  n  and  TMm  n  modes  :  equation  77 

TMq  n  and  TEM  modes  :  equation  80 

TEq  n  modes  :  equation  83. 

8.00  TERMINATION 
9.00  FORMAT  STATEMENTS 

Subroutine  PARDRV 

Subroutine  PARDRV  is  also  called  by  the  integration  routine  DGEAR 
and  it  evaluates  an  N  x  N  Jacobian  matrix  of  partial  derivatives . 

2.00  DEFINE  DIELECTRIC  AND  FIRST  DERIVATIVE  VALUES 

This  section  evaluates  FDRV  and  DREL  based  on  the  present 
radial  position  used  by  the  integration  routine  DGEAR. 

2.20  EVALUATE  PARTIAL  DERIVATIVES 

Compute  one  of  the  following  first  order  partial  differen¬ 
tial  equation  systems  for  a  given  mode  type: 

TEm  n  and  TMm  n  modes  :  equation  84 

TMq  n  and  TEM  modes  :  equation  85 

TEq  n  modes  :  equation  86. 

8.00  TERMINATION 


9.00  FORMAT  STATEMENTS 
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Function  FNCTl 


The  zero  finding  routine  ZREALl  calls  FNCTl,  which  defines  the 
functions  for  which  the  roots  of  the  TEq  n,  TMq  n  (n  >  1)  and  TEM  modes 


are  found . 
1.00 


3.00 


3.50 


A. 10 


8.00 


INITIALIZATION 

Initialize  variables  used  in  the  calling  argument  of  the 
Integration  routine  DGEAR.  Also,  set  initial  conditions 
for  one  of  the  following  mode  types: 

TMq  n  and  TEM  modes  :  equation  89 


TEq  n  ®°^es 


:  equation  91. 

CALL  DGEAR 

The  integration  routine  DGEAR  is  called. 

SET  UP  FOR  NEXT  DIELECTRIC  SECTION 

If  a  TMq  n  or  TEM  mode  type  is  chosen,  set  equation  73 
DEFINE  FUNCTION  STATEMENT  FNCTl 

According  to  which  mode  type  is  chosen,  evaluate  one  of  the 
following  equations: 

TMq  n  and  TEM  modes  :  equation  93 
TEq  n  modes  :  equation  94. 

TERMINATION 


9.00  FORMAT  STATEMENTS 


Function  FNCT2 


The  N  dimensional  minimization  routine  E04JBF  calls  FNCT2 ,  which 
defines  the  function  for  the  roots  of  the  TEm  n  and  TMpj  n  (n  >  1)  modes. 
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1.00  INITIALIZATION 

Initialize  variables  ujad  in  the  calling  argument  of  the 
integration  routine  DGEAR  and  equation  87. 

3.00  CALL  DGEAR 

The  integration  routine  DGEAR  is  called. 

3.50  SET  UP  FOR  NEXT  DIELECTRIC  SECTION 
Compute  equation  72. 

4.10  DEFINE  FUNCTION  STATEMENT  FNCT2 
Evaluate  equation  95. 

8.00  TERMINATION 
9.00  FORMAT  STATEMENTS 

Subroutine  NORMAL 

Subroutine  NORMAL  normalizes  the  field  components  of  a  given  mode. 

2.00  COMPUTE  SPLINE  COEFFICIENTS 

Compute  the  spline  interpolation  of  the  field  components 

for  one  of  the  following  mode  types: 

TE_  _  and  TM  modes  :  E  ,EX ,H  ,Ha 
nn,n  m,n  r’  $  r  $ 

TMn  and  TEM  modes  :  E  ,H 
0,n  r  ip 

TEn  modes  :  E^.H  . 

u  ,n  $  r 

3.00  TRANSVERSE  COMPONENT  INTEGRATION 


Using  the  integration  routine  DCADRE  (which  calls  Function 
CXINT) ,  equation  68  is  computed. 
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NORMALIZE  FIELD  COMPONENTS 


Normalize  the  field  components  of  the  chosen  mode  by  divid¬ 


ing  by  the  square  root  of  the  computed  integral  in  section 


3.00. 


TERMINATION 


FORMAT  STATEMENTS 


Subroutine  ORTHO 


Subroutine  ORTHO  computes  the  orthogonality  between  two  chosen 


modes . 


2.00  OUTER  FREQUENCY  LOOP 


Compute  orthogonality  at  L  multiple  frequencies. 


3.00  INNER  LOOP  TO  EVALUATE  ORTHOGONALITY  BETWEEN  TWO  MODES 


Determine  which  pair  of  modes  are  to  be  evaluated. 


4.00  COMPUTE  SPLINE  COEFFICIENTS 


Compute  the  spline  interpolation  of  the  field  components 


for  one  of  the  following  mode  types: 


TE_  „  and  TM_  _  modes  :  E  ,E , ,H  ,H 
m,n  m,n  r  4  r  $ 

TMq  n  and  TEM  modes  :  Ef ,H^ 


TE0,n  modes 


:  E**H  • 
4>  r 


5.00  TRANSVERSE  COMPONENT  INTEGRATION 


Using  the  integration  routine  DCADRE  (which  calls  Function 


CXINT)  on  equation  68,  orthogonality  is  computed. 


6.00  END  OF  INNER  LOOP 


7.00  END  OF  OUTER. LOOP 
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8.00  TERMINATION 


9  .00  FORMAT  STATEMENTS 


Function  CXINT 


Function  CXINT  is  called  by  the  integration  routine  DCADRE  and 

defines  the  argument  of  Eq.  68. 

2.00  EVALUATE  SPLINE  COEFFICIENTS 

Evaluate  the  spline  coefficients  for  the  field  components 

of  one  of  the  following  mode  types: 

TE  and  TM_  _  modes  :  E  ,E  ,H  ,H 
m,n  m  ,n  r  $  r  $ 

TMn  and  TEM  modes  :  E  ,H 
u  ,n  r  <p 

TEn  _  modes  :  E,,H  . 

U,n  r 

3.00  DEFINE  ARGUMENT  FOR  NORMALIZATION 

Based  on  whether  k_  is  real  (propagating)  or  imaginary 

z 

(nonpropagating)  ,  compute  one  of  the  following  arguments 


for  equation  68: 


k  real 
z 


:  E  H  -  EH 
r  ip  <p  r 


k2  imaginary  :  E^H.  -  E^H^  . 

A. 00  DEFINE  ARGUMENT  TO  COMPUTE  ORTHOGONALITY 

Compute  arguments  of  equation  68  based  on  one  of  the  fol¬ 
lowing  cases: 

4.10  MODE  1  PROPAGATING,  MODE  2  PROPAGATING 

k  for  both  modes  is  real;  define  EH-  EH 


r  $  <p  r 


4.20  MODE  1  PROPAGATING,  MODE  2  CUTOFF 


k2  for  mode  1  is  real  and  imaginary  for  mode  2;  define 

EaH  -EH, 

♦  r  r  $ 
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4.30  MODE  1  CUTOFF,  MODE  2  PROPAGATING 


kz  for  mode  1  is  imaginary  and  real  for  mode  2:  define 

EH.  -  E.H 
r  p  p  r 

MODE  1  AND  MODE  2  CUTOFF 

kz  is  imaginary  for  both  modes;  define  E^H^  -  E^H^ 


TERMINATION 
FORMAT  STATEMENTS 


Subroutine  DPLOT 


Subroutine  DPLOT  generates  the  field  plots  for  a  given  mode  which 
can  be  processed  into  making  hard  copy  printouts. 


4.2  Process  for  a  General  Dielectric  Profile 


The  two  major  tasks  in  the  program  concern  finding  the  propagation 
constant  kz  (and  where  appropriate,  a)  and  corresponding  field  compo¬ 
nents  for  a  desired  mode.  To  accomplish  this,  the  second  order  differ¬ 
ential  equation  system  (Eqs.  47  or  48)  must  be  integrated  in  the  radial 
coordinate  across  the  dielectric  profile.  For  a  dielectric  profile  free 
of  any  discontinuities,  there  is  a  single  integration  from  the  inner  to 
outer  conductor  radius.  If  there  are  N  discontinuities,  then  integra¬ 
tion  must  be  performed  N  +  1  times. 

The  dielectric  profile  is  originally  defined  at  discreet  points. 
As  the  integration  is  being  performed  from  the  inner  to  outer  conductor 
radius,  a  point  may  be  chosen  by  the  integration  routine  at  which  the 
dielectric  is  not  defined.  Thus,  before  the  integration  is  carried  out, 
cubic  spline  interpolation  is  performed  to  generate  a  smooth  curve  for 
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the  dielectric.  A  discussion  involving  the  theory  of  cubic  splines  used 
in  the  program  is  presented  by  de  Boor.4  The  interpolatory  subroutines 
used  are  presented  in  Appendix  A  (Section  A. 5).  Once  the  Interpolation 
has  been  completed,  integration  may  proceed  as  before.  For  a  dielectric 
profile  free  of  any  discontinuities,  the  standard  cubic  spline  interpo¬ 
lation  is  performed.  If  the  dielectric  profile  is  discontinuous,  the 
standard  spline  approach  generates  a  large  "ringing"  in  the  vicinity  of 
the  jump.  The  approach  taken  in  solving  this  problem  is  described  in 
Appendix  C. 

Once  kz  has  been  determined  for  a  given  mode,  the  field  components 
as  a  function  of  r  need  to  be  generated.  The  integration  routine  used 
to  integrate  the  field  components  from  the  inner  to  outer  conductor  is 


the  same  as  that  used  in  the  “shooting"  method  to  find  k  (and  a). 

z 

Again,  cubic  spline  interpolation  must  be  used  on  the  dielectric  profile 
in  order  to  carry  out  the  integration. 


4 . 3  Finding  a  Mode 


In  order  to  find  the  mode  propagation  constant  kz  and  resulting 


field  patterns  for  a  given  mode,  an  initial  guess  for  kz  must  be  pro¬ 
vided.  Whenever  a  dielectric  different  from  air  is  introduced  into  a 
waveguide,  the  cutoff  frequency  of  a  given  mode  is  lowered  with  respect 
to  that  of  air.  This  is  due  to  the  fact  that  electromagnetic  waves 
travel  at  a  reduced  velocity  governed  by  the  following  equation: 


0 

K 


(96) 
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where 


« 


c£  ■  speed  of  light  in  dielectric 

Cq  *  speed  of  light  in  air 

*  relative  dielectric  permittivity 

With  this  fact  in  mind,  the  initial  guess  for  kz  falls  between  the 

values  of  k  in  the  limiting  cases  of  e  =  e  =*  1  (air)  and  e  ,  *  e  „ 

z  6  rl  r2  rl  r2 

=  desired  relative  dielectric  permittivity, 


where 

k  “  mode  propagation  constant  for  e  «*  e.  *  1 
za  r  rl  r2 

kzh  “  mode  propagation  constant  for  *»  desired  dielectric 

value 

As  mentioned  in  Section  3.2,  an  additional  guess  for  the  dimensionless 

parameter  a  is  required  for  TE„  _  (HE  )  and  TM_  _  (EH  )  modes. 

m,n  v  n.n'  Tm,n 

To  evaluate  the  geometry  of  Fig.  3A  for  a  chosen  dielectric  pro¬ 
file,  one  must  "build  up"  to  the  final  desired  dielectric  profile.  This 
is  achieved  by  starting  with  e  •  ■  1  and  gradually  perturbing  the 

system  by  Increasing  the  value  of  e  j  and/or  e  .  This  is  schematically 
shown  in  Fig.  4.  The  initial  guess  for  kz  at  the  Mth  step  is  based  on 
Eq.  97,  where  kzh  18  the  mode  propagation  constant  for  the  homogeneously 
filled  coax  with  the  current  dielectric  value.  For  the  ra  *  0  modes,  the 
Initial  guess  for  a  at  the  Mth  step  is  set  equal  to  the  final  value  from 
the  previous  step,' 
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When  m  ■  0,  a  one-dimensional  zero  finding  routine  (ZREAL1,  Section  A. 2) 
is  used  to  find  the  zero.  When  m  *  0,  a  two-dimensional  minimum  finding 
routine  (E04JBF,  Section  A. 3)  is  used  to  find  the  minimum.  In  this 
second  case,  a  relatively  poor  initial  guess  for  kz  or  a  will  nudge 
E04JBF  towards  an  undesired  minimum.  In  particular,  the  initial  guess 
for  a  was  critical.  This  result  was  found  empirically  by  running  the 
program  and  observing  the  affects  of  varying  a.  Conversely,  for  the  ra  * 
0  cases,  the  initial  guess  for  kz  was  not  as  critical. 

To  have  any  confidence  that  the  final  value  of  kz  at  the  Mth  step 
represents  the  desired  mode,  certain  conditions  are  examined.  The  first 
condition  requires  that  kz  fall  between  the  limits"  of  kza  and  kzh,  as 
specified  by  Eq.  97.  If  this  is  not  the  case,  then  the  program  has 
unwittingly  ''walked”  to  another  solution.  The  next  check  involves 
examining  the  evolution  of  the  field  plots  from  the  air  filled  case  to 
the  present  dielectric  profile.  Characteristics  which  should  be  scruti¬ 
nized  include: 

1.  The  number  of  zero  crossings 

2.  Conformance  to  the  boundary  condition  that  the  tangential 
components  of  the  electric  field  E  and  the  normal  components 
of  the  magnetic  field  H  at  r  ■  a  and  r  *  c  are  zero 

3.  Relative  magnitudes  of  similar  field  components 

4.  Poynting’s  vector 

The  last  major  criterion  to  be  satisfied  is  orthogonality,  which  is 
computed  between  the  current  mode  of  interest  and  another  conveniently 
chosen  mode. 
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V.  PRESENTATION  OF  RESULTS 


The  geometry  of  Fig.  3A  was  examined  for  the  following  modes: 
TMq  j,  TEq  j,  TEM,  TMj  ^(EH^  ^).  The  dimensions,  dielectric  profile 
values,  and  frequency  used  for  each  mode  were 


a 

b 

c 

e 


rl 


'r2 


1.0  cm 
1.5  cm 
2 .0  cm 
1.0,  2.0 
1.0 


There  are  two  cases  presented  for  each  mode: 


Case  1: 


'rl 


r2 


1.0 


Case  2:  Scl  -  2.0,  -  1.0 


Both  cases  were  evaluated  at  32  gigahertz  (GHz).  For  clarity,  region  1 
will  refer  to  the  dielectric  region  of  e  ,  and  region  2  will  refer  to 
the  dielectric  region  of 

The  total  number  of  points  at  which  the  dielectric  profile  was 
defined  was  46.  The  points  were  equally  spaced,  except  in  the  neighbor¬ 
hood  of  the  discontinuity  (b  -  1.5  cm).  Here,  points  were  concentrated 
in  an  effort  to  increase  the  accuracy  of  the  integration  and  when  nor¬ 
malization  and  orthogonality  were  computed.  In  determining  the  optimum 
number  of  points  to  use,  consideration  had  to  be  given  to  the  cost, 


time,  and  numerical  accuracy  desired.  The  number  of  points  chosen,  46, 


gave  good  numerical  accuracy 

while  keeping  the  time 

and 

cost  to  a  mini-  i 

mum.  Appendix  E  lists 

the  values  of  the  individual 

points. 

Before  proceeding 

,  the 

results  for  each  mode 

are 

presented  below 

(e^2  is  always  equal  to 

!)• 

Mode 

fri 

k  (meters  M 

a 

™0,1 

1.0 

592.99 

2.0 

681.37 

cs/ 

TE0,1 

1.0 

589.09 

2.0 

819.66 

/V 

TEM 

1.0 

670.89 

2.0 

916.37 

™1,1 

1.0 

589.84 

-3.72  E-8 

EHj  j  mode  1 

2.0 

678.27 

-7.33  E-2 

EHi  ^  mode  2 

2.0 

815.32 

-4.11 

Sections  5.1  through 

5.4  will  corapare/constrast 

the 

field  component 

plots  for  the  TMq  j , 

TE0,1 

,  TEM,  and  TMj^  (E^  J 

modes ,  respec- 

tively.  Computation 

of  orthogonality  will  be  presented  in  Section 

5.5.  Two  dispersion  plots  will  then  be  presented  and  discussed  in 
Section  5.6.  Finally,  Section  5.7  will  show  the  effects  on  the  propaga¬ 
tion  constant  k  as  the  dielectric  ratio  of  e  ,  to  e  is  varied, 
z  rl  r2 
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5.1  THn  Mode 

The  effect  of  the  discontinuous  dielectric  profile  at  r  *  b  is 
readily  apparent  in  Fig.  5  for  Er.  Notice  that  the  magnitude  of  Er 
doubles  (as  required  for  an  inner  to  outer  dielectric  ratio  of  two) 
along  with  the  fact  that  the  field  lies  primarily  in  region  2.  As  with 
the  air  filled  case  (Fig.  6),  there  is  one  zero  crossing  in  each  case. 

The  Ez  component  distribution  has  evolved  from  being  almost  sym¬ 
metric  (Fig.  7)  to  primarily  concentrated  in  region  1  (Fig.  8).  Note 
that  Ez  is  continuous  at  r  *  b  (Fig.  8),  but  its  first  derivative  with 
respect  to  r  is  discontinuous.  The  boundary  condition  that  Ez  be  equal 
to  zero  at  r  =»  a  and  r  =  c  is  satisfied  in  both  plots. 

Figures  9  and  10  depict  the  H,  component  for  cases  1  and  2, 

9 

respectively.  Case  2  (Fig.  9)  shows  H  continuous  at  r  *  b,  with  the 

9 

overall  field  distribution  primarily  in  region  1  (as  with  Fig.  10). 

Poynting's  vector  shows  that  for  case  2  (Fig.  11),  the  energy 


resides  primarily  in  region  2,  implying  that  most  of  the  mode  propagates 
in  this  region.  This  boldly  contrasts  the  energy  distribution  for  case 
1  (Fig.  12),  which  is  slightly  larger  for  region  1. 
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Flf.  7.  THq  j  •o4«  at  32  Cfla,  MI  •  1;  ■  coaponant. 
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Fig.  8.  aoda  at  32  CHa,  EM  »  1;  Pointing’ a  vactor  ErH^. 
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Fig.  9..  TMq # 1  aod«  *t  32  CHx ,  ESI  «  2;  lr  component. 


'  ft 

A 


POYNTINC’S  VECTOR-  ERwHT 


9.UQC-C3 


1.S0C-C2  1 

IVDM.  rOSZTION  -f 


a.ioc-oa 


2.40C-42 


Fig .  11.  TMqj  ao4«  at  32  CHs ,  Ell  •  2;  component. 
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5.2  TEq  x  Mode 

The  E.  component  shows  that  the  field  has  shifted  primarily  to 
9 

region  1  (Fig.  13)  from  what  was  almost  a  symmetrical  distribution  (Fig. 

14).  In  both  cases,  E  is  equal  to  zero  at  r  ■  a  and  r  »  c,  and  is 

9 

continuous  at  r  ■  b. 

The  Hr  component  is  similarly  affected  in  that  most  of  the  field 
lies  in  region  1  (Fig.  15)  in  contrast  to  the  air  filled  case  (Fig. 
16) .  The  boundary  condition  that  Hr  be  equal  to  zero  at  r  -  a  and  r  *  c 
is  satisfied  in  both  plots. 

The  Hz  component  for  case  2  (Fig.  17)  continues  the  trend  toward 
concentrating  in  region  1.  Although  continuous  at  the  dielectric  dis¬ 
continuity,  the  first  derivative  of  H2  with  respect  to  r  is  discontinu¬ 
ous.  Also ,  there  is  one  zero  crossing,  as  with  the  air  filled  case 
(Fig.  18). 

Poynting's  vector  for  case  2  (Fig.  19)  shows  that  the  majority  of 
the  energy  is  in  region  1  compared  with  the  almost  symmetrical  distribu¬ 
tion  for  case  1  (Fig.  20).  Consequently,  the  TEq  j  mode  for  the  dielec¬ 
tric  profile  of  case  2  rropagates  primarily  in  region  1. 
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Fig.  13.  TEG,!  aoda  »t  32  CHx,  ER1  ■  1;  I  coaponant. 
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5.3  TEM  Mode 


Comparing  Figs.  21  (case  1)  and  22  (case  2)  for  the  Er  component 
reveals  that  a  greater  percentage  of  the  field  Ilea  In  region  1  for  case 
2.  The  appropriate  discontinuous  jump  at  r  -  b  la  satisfied  for  case  2. 
Theoretically,  the  TEM  mode  has  no  E  (or  H  J  component  when 

T. 


t  ,  •  i  ,  -  1.  Figure  23  shows  values  for  the  air  filled  case  which, 

r  1  r2 


although  not  exactly  iero,  are  very  email  and  reflect  the  accuracy  of 
the  program.  But,  the  values  for  Fig.  24  show  unquestioned  existence  of 
the  Ej,  component  when  t  -  2.  Here,  a  peak  Is  reached  at  the  dielec¬ 
tric  Interface  (r  «  b)  where  the  field  Is  continuous,  but  Its  first 


derivative  with  respect  to  r  is  discontinuous.  As  required  E^  Is  zero 


at  r  *  a  and  r  ■  c. 


The  H  component  shows  an  increased  concent  at  1  on  In  region  1 


'hlg.  25)  relative  to  that  of  the  air  filled  case  (Fig.  26).  At  r  - 


b.  H  Is  continuous,  although  Its  first  derivative  with  respect  to  r  Is 


d 1  sc  on  t 1 nuous . 

Poyntlng’s  vector  shows  how  the  mode  has  evolved  Into  propagating 
almost  entirely  In  region  1  (Fig.  27)  compared  to  that  of  the  air  filled 
case  (Fig.  26). 
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Fig.  21.  TIM  aod«  *t  32  GHz  ,  Ell  ■  l;  E  component. 
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Fig.  2*.  TEM  aod«  it  32  GHt ,  EH  -  1;  Poyntlng'a  vector  ErH^. 
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Pl|.  25.  TEM  «od*  at  32  GH* ,  ER1  -  2;  Er  component. 
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Fig.  27.  TEM  Mode  «C  32  CH* ,  ER1  *  2;  component. 
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Ftg.  28.  TEM  sod*  at  32  CHt ,  ERI  »  2;  Poyntlng’e  vector  EfH^. 


5.4  TM  vEHj  A)  Mode 

Here,  there  are  two  sets  of  answers  when  c  ,  -  2  signlrving  a 

r  1 

split  of  the  TM,  .  mode  upon  hybridization.  The  plots  for  the  hvbru 
ized  modes  are  labeled  EH^  ^  mode  1  and  EH^  ^  mode  2  correspond  ir.g  to 
the  values  of  b78.27  and  815.32.  respectively.  In  the  discussion,  we 
will  simply  refer  to  them  as  mode  1  and  mode  2. 

The  Ef  component  for  modes  1  and  2  iFigs.  29  and  3  ' ,  respectively' 
undergoes  the  required  jump  at  r  ■  b  by  doubling  in  magnitude.  Both 

modes  have  one  zero  crossing,  as  with  the  TM^  ^  mode  (Fig.  31'.  But, 
mode  1  has  its  field  distribution  concentrated  in  region  2,  while  mode  2 
has  its  field  almost  entirely  in  region  1.  Note  that  the  Er  component 
for  mode  1  is  from  one  to  two  orders  of  magnitude  larger  than  mode  2 

over  most  of  the  radial  cross  section. 

The  azimuthal  electric  field,  E  ,  for  mode  1  (Fig.  32)  has  evolved 

♦ 

into  two  positive  peaks  compared  to  the  one  positive  peak  for  mode  2 
(Fig.  33).  Both  modes  1  and  2  have  their  fields  concentrated  in  region 

1,  as  with  the  TM.  .  mode  (Fig.  34).  In  all  the  plots,  E  is  continuous 

i  ,  i  <9 

at  r  -  b  and  is  zero  at  r  *  a  and  r  -  c.  We  further  note  that  the  t 

i 

component  for  mode  2  is  one  to  two  orders  of  magnitude  larger  than  model 
over  the  radial  cross  section. 
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The  E_  component  shows  that  for  both  modes  1  and  2  (Figs.  35  and 
36,  respectively),  the  field  resides  primarily  in  region  1.  Note  that 
mode  1  goes  through  a  zero  crossing.  Along  with  the  TMj  ^  mode  (Fig. 
37),  all  three  plots  satisfy  the  boundary  condition  that  Ez  be  equal  to 
zero  at  r  ■  a  and  r  -  c  as  well  as  being  continuous  at  r  «*  b.  We  again 
point  out  that  mode  1  is  an  order  of  magnitude  larger  than  mode  2  over 
the  radial  cross  section.  This  is  consistent  with  the  fact  that  the 
magnitude  of  a  (Eq.  49)  for  mode  1  is  less  than  that  for  mode  2,  result¬ 
ing  in  Ez  being  larger  for  mode  1. 

Figures  38,  39,  and  40  show  the  component  for  the  TMj  ^  mode, 
mode  1  and  mode  2,  respectively.  We  note  the  field  distribution  having 
an  increased  concentration  in  region  1  for  modes  1  -(which  goes  through  a 
zero)  and  2.  The  Hr  component  for  mode  2  is  as  much  as  two  orders  of 
magnitude  larger  than  mode  1  over  the  radial  cross  section.  At  r  -  a 
and  r  «  c,  Hr  is  equal  to  zero  for  all  three  modes,  as  required. 

As  with  the  H_  component,  H,  shows  a  majority  of  the  field  resid- 
r  $ 

ing  in  region  1  for  modes  1  and  2  (Figs.  41  and  42,  respectively).  Note 

that  mode  2  does  not  have  a  zero  crossing  unlike  mode  1  and  the  TMj  j 

mode  (Fig.  43).  All  three  modes  are  continuous  at  r  «  b.  Here,  the  H, 

9 

component  for  mode  1  is  an  order  of  magnitude  larger  than  mode  2  over 
the  radial  cross  section. 
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The  fact  Chat  Che  magnitude  of  a  is  larger  for  mode  2  than  for 
mode  1  is  evidenced  in  the  Hz  component  (Figs.  44  and  45.  respectively). 
The  field  distributions  are  both  concentrated  in  region  1,  but  the  mag¬ 
nitude  of  H  for  mode  2  is  larger.  Along  with  the  TM,  .  mode  (Fig.  46), 

all  three  modes  have  one  zero  crossing  and  are  continuous  at  r  *  b. 

Finally,  the  overall  affects  of  the  dielectric  profile  on  modes  1 
and  2  are  summarized  by  examining  the  cosine  (e^H^)  and  sine  (E^H^) 
terms  of  Poynting's  vector.  The  cosine  term  for  the  TM^  ^  mode  (Fig. 
47)  shows  initially  that  most  of  the  energy  lies  in  region  1.  The  sine 
term  for  the  TMj  j  mode  (Fig.  48)  shows  a  slightly  greater  concentration 
of  energy  in  region  1.  Similar  to  the  TM^  1  mode,  the  cosine  term  for 
mode  2  (Fig.  49)  has  an  even  larger  fraction  of  energy  residing  in 
region  1,  and  its  sine  term  (Fig.  50)  has  almost  all  of  its  energy  in 
region  1.  However,  mode  1  shows  behavior  that  contrasts  that  of  the 
TMj  ^  mode  and  mode  2.  Figure  51  for  mode  1  shows  that  the  majority  of 

the  cosine  energy  term  resides  in  region  2.  The  sine  term  of  mode  1 

(Fig.  52)  has  its  energy  primarily  in  the  peak  of  region  1,  but  the 
magnitude  is  negative. 
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Fig.  50.  node  2  ac  32  GHz,  ERI  ■  2;  Hz  component. 
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From  examining  the  field  component  plots  and  the  sine  and  cosine 
Poynting  vector  components  (which  combine  to  represent  the  total  net 
flow  of  energy),  we  note  that  the  TM^  ^  mode  and  mode  2  have  similar 
characteristics.  While  mode  1  is  also  similar  to  the  TM^  ^  mode  in  some 
respects,  its  energy  distribution  is  quite  different.  Since  sine  and 
cosine  are  functions  which  are  90  degrees  out  of  phase,  there  will  be 
(for  a  given  cross-sectional  plane)  points  in  the  azimuth  where  cosine 
is  zero  and  sine  is  not.  Although  the  majority  of  energy  for  mode  1 
propagates  in  region  2  in  the  +z  direction,  there  are  locations  in  the 
cross-sectional  plane  where  there  is  a  net  flow  of  energy  (~1.0  E-4)  in 
the  -z  direction  (as  a  result  of  the  sine-cosine  relationship).  This 
behavior  is  analogous  to  the  "back  eddys"  created  from  water  flowing 
through  a  pipe  containing  obstacles.  Unlike  mode  1,  the  energy  for  the 
TM1  j  mode  and  mode  2  propagates  primarily  in  region  1  in  the  +z  direc¬ 
tion. 


5 . 5  Evaluation  of  Orthogonality 

Sections  5.1  through  5.4  presented  the  values  for  kz,  a  (where 
appropriate),  and  field  plots  for  the  given  modes.  With  every  mode,  kz 
for  the  dielectric  profile  of  e  ■  2  (e^j  -  l)  was  larger  than  that  for 
the  air  filled  case.  This  is  consistent  with  Eq.  97  of  Section  4.2. 
From  Section  5.4,  the  magnitude  of  o  showed  that  significant  hybridiza¬ 
tion  had  occurred  for  modes  1  and  2.  The  fact  that  a  for  the  TM^  j  mode 
was  finite  (but  very  small)  reflects  the  numerical  accuracy  of  the 
program.  The  boundary  conditions  at  r  »  a  and  r  »  c  for  the  tangential 


electric  field  and  normal  magnetic  field  components  were  satisfied  for 


all  the  modes  examined.  The  boundary  condition  that  Er  jump  by  a  factor 

of  two  f ratio  of  e  ,  to  e  .1  at  r  ■  b  for  case  2  was  also  satisfied. 
v  rl  r2' 

Poynting’s  vector  has  shown  for  each  mode  that  a  discontinuous  dielec¬ 
tric  profile  shifts  the  majority  of  energy  to  one  of  the  two  regions. 
Thus,  as  previously  noted,  the  mode  propagates  primarily  in  one 
region.  One  further  aspect  that  must  be  examined  is  whether  the  modes 
in  question  are  orthogonal. 

Orthogonality  was  computed  between  the  following  pairs  of  modes: 

TEM  and  TMq  TEq  ^  and  TEq  j  (chosen  for  convenience);  EH^  ^  mode  1, 

and  EHj  j  mode  2.  The  results  of  Eq.  68  (Section  3-5)  for  each  pair  are 

presented  under  columns  2  through  5  where  the  mode  before  the  backlash 

contributes  the  E_  and  E,  components  and  the  mode  following  the  backlash 
r  $ 

contributes  the  H,  and  H_  components: 
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Mode  and  EH 
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Crl  Mode  1/Mode  1  Mode  1/Mode  2  Mode  2/Mode  1  Mode  2/Mode  2 

2.0  1.0  2.87  E-6  1.94  E-6  1.0 

Theoretically,  if  two  modes  are  orthogonal,  then  Eq.  68  is  exactly  equal 
to  zero.  The  results  show  that  for  each  mode  pair,  the  computed  values 
of  orthogonality  are  very  small  for  both  values  of  e  .  Furthermore, 
columns  2  and  5  (which  represent  computation  of  self  orthogonality)  have 
the  value  of  1  for  each  case  of  e  ^  ,  as  expected.  In  light  of  these 
results,  we  conclude  that  the  modes  TMq  j,  TEq  j,  TEM,  EHj  j  mode  1,  and 
EH ^  j  mode  2  are  valid. 

5.6  Dispersion  Plots 

Two  dispersion  plots,  w  versus  k  ,  are  presented  in  this  section 

z 

for  the  following  conditions: 


Frequency  range  :  1.0-34.0  GHz 

Figure  53  is  a  plot  of  the  TMQ  p  TEq ^ 1 ,  and  TEM  modes,  and  Fig.  54  is  a 

plot  of  the  modes  EH^  j  mode  1  and  EHj  ^  mode  2.  In  each  plot,  the  y 

axis  (to)  has  been  normalized  by  dividing  by  the  speed  of  light  in 

vacuum  (c^)  and  multiplying  by  the  radius  of  the  dielectric  interface,  b 

(b  ■  1.5  cm).  Upon  multiplying  the  x  axis  (k  )  by  b  results  in  both 

z 

axes  being  dimensionless  and  of  the  same  order  of  magnitude.  Thus,  the 
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x  and  y  axes  are  labeled  kzb  and  kb  (where  k  -  wc~*)»  respectively.  A 
negative  value  of  kzb  represents  the  propagation  constant  when  it  is 
cutoff  (imaginary). 

Beyond  a  certain  frequency  for  a  given  mode,  both  plots  reveal 

propagation  between  the  two  velocity  of  light  curves  v  «  cn  and  v  *  c 

u  e 

(for  the  dielectric  profile  of  -  e -  2).  For  the  mode  in  ques¬ 
tion,  this  implies  that  the  fraction  of  the  total  energy  which  propa¬ 
gates  in  region  2  is  a  slow  wave  mode  (below  v  «  c^)  and  the  remaining 
energy  which  propagates  in  region  1  is  a  fast  wave  mode  (above  v  -  C£). 

Thus,  the  mode  propagates  with  an  overall  k  that  lies  between  the  k 

z  z 

for  the  air  filled  region  (region  2,  v  «  c  )  and  the  k  for  region 

0  z 

2  (erj  “  2 ,  v  «  c^).  This  is  precisely  Eq.  97  presented  in  Section 
4.2.  Hence,  the  dispersion  plots  reinforce  the  validity  of  the  modes 
examined  for  the  dielectric  profile  of  -  2  and  *  1. 

A  connection  can  be  made  between  the  asymmetric  behavior  of  the 
mode  as  the  frequency  goes  to  infinity  with  the  distribution  of  the 
energy  from  the  analysis  of  Poynting’s  vector.  For  a  mode  that  is 
concentrated  in  the  dielectric  (tEq  j,  TEM,  EH^  ^  mode  2),  the  disper¬ 
sion  should  approach  the  v  ■  c£  line.  Conversely,  for  a  mode  that  is 
concentrated  in  the  vacuum  region  (tMq  ^ ,  EHj  j  mode  l) ,  the  dispersion 
should  approach  the  v  •*  Cq  line.  This  appears  to  be  the  trend  for  the 
data  in  Figs.  53  and  54. 
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5-7  Effect  of  the  Dielectric  Profile  on  k 


The  effect  on  the  propagation  constant  k^  by  varying  the  dielec¬ 
tric  ratio  of  e  to  t  from  one  to  two  at  32  gigahertz,  is  examined  in 
this  section.  As  in  Figs.  53  and  54,  kz  is  normalized  by  multiplying  by 
b. 

Figure  55  shows  how  modes  1  and  2  split  as  the  dielectric  ratio 
increases.  Although  the  modes  are  orthogonal  to  one  another  at  each 
value  of  the  dielectric  ratio,  there  is  generally  a  slow  increase  in  the 
magnitude  of  the  computed  values  as  the  ratio  becomes  larger. 
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This  degradation  is  principally 

due  to  the  "ringing" 

which  becomes 

larger  at 

the  dielectric 

discontinuity. 

The  consequence 

is  an  increase 

in  the  errors  of  the  solutions  which,  in 

turn,  affects  the 

orthogonality 

computations . 
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Fig.  55.  K_  versus  ER1/ER2 


VI.  CONCLUSION 


6.1  Advantages  and  Drawbacks 


The  direct  integration  approach  formulated  in  Chapter  3  has  vari¬ 
ous  advantages  over  the  traditional  modal  expansion.  First,  the  TEM 
mode  can  be  examined  for  a  discontinuous  dielectric  profile  (as  was  done 
for  Fig.  3A) .  The  second  order  differential  equation  systems  for  the 
m  *  0  and  m  ■  0  modes  (Eqs.  47  and  48,  respectively)  were  derived  inde¬ 
pendently  of  the  dielectric  profile.  As  a  result,  each  differential 
equation  system  can  evaluate  a  general  dielectric  profile  (which  could 
be  a  combination  of  stepped  and  linearly  graded  sections).  This  charac¬ 
teristic  is  quite  unlike  the  modal  expansion  approach  where  the  size  of 
the  dispersion  determinant  varies  as  the  square  of  the  number  of  steps , 
and  therefore  the  determinant  grows  in  complexity  and  size  in  accordance 
with  the  complexity  of  the  profile.  Unlike  the  dispersion  determinant, 
Eqs.  47  and  48  involve  no  Bessel  functions.  Hence,  the  computer  program 
does  not  have  to  evaluate  Bessel  function  expansions,  therefore  allowing 
the  investigation  of  modes  that  are  near  cutoff  (nonpropagating).  The 
success  of  this  technique  is  illustrated  by  the  solutions  obtained  in 
Chapter  5. 

As  noted  in  Section  4.3,  the  m  *  0  modes  (TE„  _  and  TM  modes) 

’  '  m  ,n  m,n  ' 

had  a  tendency  to  "walk"  to  an  undesired  solution  for  kz  and  a  if  the 
initial  guesses  were  relataively  poor.  Also,  the  time  involved  in 
finding  a  solution  to  an  m  *  0  mode  was  a  factor  of  five  or  more  slower 
than  the  m  -  0  modes. 


£ 
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6.2  Future  Goals  and  Applications 

The  formulation  in  Chapter  3  was  derived  with  the  permeability  of 
the  dielectric  to  be  that  of  air  (Wq)-  A  natural  extension  would  be  to 
incorporate  a  permeability  profile  y(r).  Also  of  interest  would  be  the 
addition  of  a  complex  permittivity  e  *  ( r)  to  the  formulation  and  inves¬ 
tigating  the  complex  propagation  constant  kz  for  various  lossy  dielec¬ 
tric  profiles.  Finally,  work  needs  to  be  done  on  the  zero  finding  logic 
to  improve  the  speed  of  convergence  of  kz  for  a  for  the  m  *  0  modes. 
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APPENDIX  A 


NUMERICAL  SUBROUTINES 

The  following  subroutines  were  provided  by  the  IMSL5  and  NAG6 
libraries,  which  are  collections  of  mathematical  and  statistical  sub¬ 
routines  written  in  FORTRAN.  Following  the  name  of  each  routine  will  be 
the  library  from  which  it  was  taken. 

The  routine  DGEAR  is  used  to  integrate  the  second  order  differen¬ 
tial  equation  systems  (Eqs.  47  and  48).  The  zero  finding  in  the  shoot¬ 
ing  method  is  performed  by  ZREAL1  and  E04JBF.  Evaluating  the  integral 
(Eq.  68)  of  the  solutions  to  perform  normalization  and  the  orthogonality 
tests  is  carried  out  by  DCADRE.  The  routines  ICSEVU  and  DCSEVU  are  u6ed 
as  utilities  to  generate  and  evaluate  spline  approximations  to  the  given 
dielectric  profile  and  the  final  solutions  to  the  field  components. 

A.l  DGEAR  (IMSL) 

This  routine  is  used  to  integrate  Eqs.  47  and  48  in  the  "shooting 
method."  The  solution  to  a  system  of  first  order  ordinary  differential 
equations  of  the  form  y'  *  f(x,  y)  with  initial  conditions  can  be  solved 
by  DGEAR.  The  basic  methods  taken  in  obtaining  solutions  are  of  the 
implicit  linear  multistep  type.  There  are  two  classes  of  such  methods 
available  to  the  user.  The  first  is  the  implicit  Adam's  methods  (up  to 
order  twelve) ,  and  the  second  is  the  backward  differentiation  formula 
methods  (up  to  order  five)  also  known  as  Gear's  stiff  methods.  We  used 
the  second  method.  With  either  case,  an  algebraic  system  of  equations 
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oust  be  solved  at  each  step  where  a  variety  of  corrector  iteration 
methods  are  available  for  use. 

To  evaluate  the  first  order  differential  equation  system,  DGEAR 
must  call  on  two  subroutines  provided  by  the  user.  The  first  sub¬ 
routine,  DER1V,  defines  and  evaluates  the  first  order  differential 
equation  system  Y^.Y^,  Y^  given  N  (the  number  of  first  order  dif¬ 

ferential  equations),  x  (the  values  at  which  to  evalute  the  equations, 
and  Yj,  Y2,  •••,  Yn  (the  integration  variables).  The  second  subroutine, 
PARDRV,  defines  and  evaluates  the  Jacobian  matrix  of  partial  deriva¬ 
tions. 

A. 2  ZREALl  ( IMSL) 

The  routine  ZREALl  is  used  to  find  the  zero  for  the  TEM,  TMq  n, 
and  TEq  n  type  modes.  This  routine  finds  the  N  real  zeros  of  a  single 
argument,  real  function  subprogram  F(X)  which  is  supplied  by  the  user. 

Upon  supplying  X  with  N  initial  guesses  Xj ,  X2 . Xn,  the  subroutine 

uses  Muller' 6  method  to  locate  the  N  real  zeros  of  F(X) .  The  solutions 
to  F(X)  ■  0  are  returned  in  X. 

A. 3  EOAJBF  (NAG) 

The  routine  EOAJBF  is  used  to  find  the  minimums  (k  and  a)  for 

v  z  ' 

TMjjj^  ( EH^  n)  and  TEm>n  (HE^  n)  type  modes.  This  routine  employs  a 
comprehensive  quasi-Newton  algorithm  for  finding: 

1.  An  unconstrained  minimum  of  a  function  of  several  variables. 

2.  A  minimum  of  a  function  of  several  variables  subject  to  fixed 
upper  and/or  lower  bounds  on  the  variables.  No  derivatives 


are  required,  but  the  user  may  specify  continuous  first  and 
second  derivatives  (the  routine  will  usually  work,  when  there 
are  occasional  discontinuities). 


A  function  of  N  variables  F  fx,  ,  x- ,  x  1  is  minimized  subject 

12  n 

to  the  constraint, 


LJ  ‘  XJ  ‘  U3 


for  j  *  1,  2,  ...,  N,  where  Lj  is  the  lower  bound  and  Uj  the  upper 
bound.  The  user  must  specify  a  starting  point  and  an  external  function 
subroutine  FUNCT  to  calcualte  the  value  of  F(X)  at  any  point  X  in  N 
dimensional  space,  where  X  *  (x^ ,  x^ ,  ...,  xn) •  The  function  subroutine 
FUNct  defines  and  evaluates  the  function  which  is  minimized. 


FC  -  F(XC(1),  XC(2) ,  ...,  XC(n) ) 


where  XC  is  an  array  of  dimension  N,  which  contains  the  current  point  of 
evaluation.  Special  variables  that  need  to  be  defined  are: 

STEPMAX  —  specifies  an  estimate  of  the  Euclidean  distance  between 
the  solution  and  starting  point 

FEST  —  specifies  an  estimate  of  the  function  value  at  the 

minimum 

IBOUND  —  specifies  whether  the  problem  was  constrained  or 
bounded 

BL  —  array  of  dimension  N  containing  the  fixed  lower  bounds 

Lj 
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BU  —  array  of  dimension  N  containing  the  fixed  lower  bounds 

UJ 

The  user  must  also  supply  the  subroutine  MONIT  with  proper  parameter 
list.  If  desired,  MONIT  can  be  used  to  monitor  the  minimization 
process.  Subroutine  E04HBF  (from  the  NAG  library),  which  computes  the 
finite  difference  intervals  for  input  to  E04JBF,  was  modified.  Nor¬ 
mally,  using  machine  accuracy  in  its  computations,  the  accuracy  of 
E04HBF  was  altered  to  that  used  in  subroutine  DGEAR. 

A. 4  DC AD RE  (IMSL) 

The  routine  DCADRE  is  used  to  integrate  the  solutions  to  Eqs.  47 
and  48  to  normalize  the  field  components  and  to . check  orthogonality. 
Numerical  integration  of  a  function  using  cautious  adaptaive  Romberg 
extrapolation  is  performed.  In  many  instances,  DCADRE  can  handle  jump 
discontinuities.  The  user  must  supply  a  single  argument,  real  function 
subprogram  F(X). 

A. 5  ICSCCU,  ICSEVU ,  DCSEVU  (IMSL) 

The  routines  ICSCCU,  ICSEVU,  and  DCSEVU  are  presented  together, 
since  they  are  all  involved  in  the  cubic  spline  interpolation  of  a  given 
set  of  points.  The  interpolatory  approximation  to  a  set  of  points  by  a 
cubic  spline  is  performed  by  ICSCCU.  The  endpoint  conditions  are  deter¬ 
mined  automatically.  Input  to  the  routine  requires  the  number  of  points 
N,  a  set  of  points  x^  ,  where  <  xi+1  for  i  -  1,  2,  ...,  N,  and  a 
corresponding  set  of  y^  (functional)  values  for  i  »  1,  2,  ...,  N. 


Evaluation  of  the  spline  coefficients  generated  by  ICSCCU  is  performed 


by  ICSEVU.  Input  to  this  routine  requires  the  set  of  spline  coeffici¬ 
ents  and  points  where  the  spline  coefficients  are  to  be  evaluated. 
Evaluation  of  the  first  and/or  second  derivatives  of  a  cubic  spline  is 
performed  by  DCSEVU.  Input  to  this  routine  requires  the  interpolated 
spline  coefficients  and  the  points  at  which  the  first  and/or  second 
derivative  should  be  evaluated. 


w_ . 


DRELT(l) 

XPAT(l) 

LTAPER 


FALSE 


Value  of  relative  dielectric  permittivity 
at  ith  radial  position 

Radial  position  of  ith  point  (meters) 

Generate  a  linearly  tapered  dielectric 
profile 


NLMOD  —  mode  control 


MC1MX 

1 

Number  of  modes  to  be 

evaluated 

MC2( J) 

1 

Azimuthal  eigenvalue 
modes 

for  TE  n 
m,n 

LTM(J) 

TRUE 

compute  k_  and  fields  for  TMn 
z  u, 

modes 

LTE(J) 

FALSE 

Compute  kz  and  fields 

for  ^O.n 

LMX(J) 

FALSE 

Compute  k  and  fields 
.  z 

modes 

for  TE_  „ 
m  ,n 

RKZG(J) 

— 

Initial  guess  for  kz 

RATIOG(J) 

— 

Initial  guess  for  o 

NLPLT  —  printing  and  plotting  control 


LPPLOT 

FALSE 

Generate  field 
C0AD7P 

plots  for  the  print  file 

LPRINT 

TRUE 

Generate  output 

data  for  COAD7P 

LDPLOT 

FALSE 

Generate  field  ; 

plots  in  subroutine  DPLOT 

LDFL 

FALSE 

ith  field  component  output  plot  generated 

LDER 

FALSE 

Field  component 

i;  Er 

LDET 

FALSE 

Field  component 

2;  E* 

LDEZ 

FALSE 

Field  component 

3=  Ez 

LDDR 

FALSE 

Field  component 

Dr 

LDDT 

FALSE 

Field  component 

5;  D* 

LDDZ 

FALSE 

Field  component 

6;  Dz 
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FALSE  Field  component  7 ;  Bf 

FALSE  Field  component  8;  B, 

9 

FALSE  Field  component  9;  Bz 


FALSE 


FALSE  Field  component  10;  Hr 

FALSE  Field  component  11;  H 

9 

FALSE  Field  component  12;  Hz 


FALSE 


FALSE 


FALSE 


Generate  plot  of  the  Poynting  vector 

ft-  T?  U 


term  E  IT 
r  4> 


Generate  plot  of  the  Poynting  vector 

4-  t?  U 


term  E,H 
<f>  r 


NLDGR  —  DGEAR  control  variables 


INDXT 


MITERT 


1.0  E-6  Next  step  size  in  x  "(independent  vari¬ 
able) 

1.0  E-8  Relative  error  bound 

1  Indicates  the  type  of  call  to  the  sub¬ 

routines  called  by  DGEAR 

1  Iteration  method  indicator 


NLZRE  —  ZREAL1  control  variables 


ITMAX 


1.0  E-5  Convergence  criterion:  a  root,  X(l),  is 
acceptable  if  ABS  (F(X(D))  <  EPS 

2  Convergence  criterion:  a  root  is  accepted 

if  two  successive  approximations  to  a 
given  root  agree  in  the  first  NSIG  digits 

100  Maximum  number  of  iterations 


NLE04  —  E04JBF  control  variables 


MAXCAL 


Number  of  function  iterations 


1.0  E-5  Accuracy  to  which  the  solution  is  desired 
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0.0 

Estimate  of  the 
minimum 

function 

value  at  the 

STEPMX 

1.0  E+5 

Estimate  of  the  Euclidean  distance  be¬ 
tween  the  solution  and  the  starting  point 

BL( 1 ,J) 

-1000 

Fixed  lower  bound 

for  k2 

BL(2,J) 

-100 

Fixed  lower  bound 

for  a 

BU(1,J) 

1000 

Fixed  upper  bound 

for  k 

z 

BU(2 , J) 

100 

Fixed  upper  bound 

for  a 

iLCAD  —  DCADRE  control  variables 

AERR  1.0  E-5  Absolute  error  of  the  zero 

ERROR  1.0  E-5  Estimated  bound  on  the  absolute  error  of 

the  zero 
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APPENDIX  C 


INTERPOLATION  OF  A  DISCONTINUOUS  FUNCTION 

As  mentioned  in  Section  4.2,  a  discontinuous  dielectric  profile 
will  generate  a  "ringing"  in  the  neighborhood  of  the  jump  if  the  stan¬ 
dard  spline  approach  is  used.  The  approach  described  below  is  a  modifi¬ 
cation  for  a  discontinuous  function. 

Let  the  range  of  the  function  be  [a,  b] ,  with  the  given  values  at 
1  <  i  <  N,  and  Xj  ■  a  and  x^  -  b.  Let  the  values  of  the  function  at 
X£  be  equal  to  Y^.  The  standard  spline  coefficients  for  the  interval 
[x^,  x  .)  are  defined  as  a^,  b^,  c^,  and  d^,  where  the  interpolated 
function  for  x  in  the  interval  is  given  by 

2  3 

Y  ■  a^  +  b^u  +  c^u  +  d^u  (C.l) 

where  u  *  x  -  x^.  Since  the  splines  normally  are  continuous  across  the 
interval  [x  x^,  [x^  x1+1),  the  following  must  hold: 

2  3 

ai-l  +  bi-l^xi  ”  Xi-l^  +  ci-l^xi  '  xi-i^  +  di-l(xi  *  xi-i^  *  ai  (C*2) 

Therefore,  a  discontinuity  may  be  introduced  across  an  interval  junction 
by  modifying  the  a^  term  only. 

Introduce  a  discontinuity  across  the  interval  junction  by  defining 
Y^  as  the  value  at  the  open  end  of  the  left  interval  and  T^  as  the  value 
at  the  closed  end  of  the  right  interval.  Also,  only  let  Tj  be  defined 
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for  i  contained  in  I,  the  set  of  interval  Junctions  that  are  discontinu¬ 
ous.  Next,  define  a  new  set  of  function  values  by 


Zi  “  Yi  ‘  l  ^  ‘  YJ  (C,3) 

J  1  2  3 


By  construction,  this  set  of  points  is  "smooth,"  since  the  discontinui¬ 
ties  have  all  been  subtracted  off.  The  standard  spline  coefficients  are 
then  calculated  for  these  values.  Finally,  the  discontinuity  is  rein¬ 
troduced  by  modifying  the  constant  terms  a^  to  A^, 

Ai  •  *i  -  ^  -  V  'c-4> 

Using  this  modified  set  of  spline  coefficients  (a^  b^  c^ ,  d^)  will 
generate  the  correct  step  discontinuity  across  the  interval  junctions 
and  a  smooth  value  for  the  function  inside  the  intervals. 


APPENDIX  D 


VARIABLE  DEFINITIONS 


The  variables  used  in  the  main  program  and  subroutines  are  pre¬ 
sented.  For  compactness,  variables  which  are  used  throughout  the  pro¬ 
gram  will  be  defined  only  in  the  section  where  they  first  appear.  The 
variables  listed  for  the  subroutines  are  local.  Section  D.7  will  define 
the  variables  set  by  the  program  and  used  in  the  calling  argument  of  the 
IMSL  and  NAG  routines.  Before  proceeding,  we  define  the  index  variables 
used  in  the  arrays. 


—  radial  position 


—  mode 


—  frequency 


D.l  Main  Program  (C0AD7R) 


Titles  and  Headers  —  used  for  output  file  purposes 


Variable 

Name 

CFTIT 

CHESS 

CVERS 


Function 


Type 


Units 


Title  for  hard  copy  field  plots  Character 

Label  for  hard  copy  field  plots  Character 

Title  used  in  the  output  print  Character 
file  C0AD7P 

Current  date  for  in  C0AD7P  - 

Current  time  of  printing  for  in 
C0AD7P 
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Variable 

Name 

Function 

Type 

Units 

Constants  — 

physical  and  numerical  constants 

P102 

»/2 

Real 

— 

EPSO 

Eg  (permittivity  of  air) 

Real 

Farad  meter-* 

UO 

Ug  (permeability  of  air) 

Real 

Henry  meter-* 

VLIGHT 

Cq  (speed  of  light  in  vacuum) 

Real 

Meters  s-* 

Program  Flow 

Control 

K 

Mode  index 

Integer 

— 

L 

Frequency  index 

Integer 

— 

FREQ 

Current  value  of  the  frequency 

Real 

GHz 

DELFRQ 

Frequency  step 

Real 

GHz 

SRKZ(K) 

Sign  of  the  final  value  of  kz 
for  the  kth  mode 

Real 

— 

Dielectric  Profile  —  The  following  variables  and  arrays 

define  and 

aid  the  processing  of  the 

dielectric  profile 

N JUMPS 

Total  number  of  discontinuities 
in  the  dielectric  profile 

Integer 

— 

NPTS 

Actual  number  of  points 
(NPTST-NJUMPS)  at  which  the 
dielectric  profile  is  defined 

Integer 

XPA(l) 

Radial  positions  at  which 
dielectric  profile  is  defined; 
length  that  of  XPAT( 1)-NJUMPS 

Real 

— 

DRELA(l) 

Defines  dielectric  profile, 
length  is  that  of  DRELT-NJUMPS 

Real 

— 

DRELN(l) 

DRELN(l)  -  LOGe  (DRELN(l)) 

Real 

— 

DRELJ(l) 

Value  of  the  "Jump”  at  the  ith 
discontinuity  in  the  dielectric 

Real 
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Variable 


Name 

Function 

Type 

DRATIO(l) 

Ratio  of  DRELA  (1+1)  to 

DRELA  (1)  at  the  ith  discon¬ 
tinuity  in  the  dielectric 

Real 

AJUMP(l) 

Radial  position  of  the  ith 
discontinuity  in  the  dielectric 

Real 

IJUMP(l) 

Index  of  the  radial  position  of 
ith  discontinuity  in  the 
dielectric 

Real 

dreldCD 

“Smoothed”  dielectric  profile; 
DRELD(l)  -  DRELT(l)  -  DRELJ(l) 

Real 

DRELDS(1 ,3) 

Spline  coefficients  for  DRELDS 
at  ith  radial  position 

Real 

DRELNS(1 ,3) 

Spline  coefficients  for  DRELNS 
at  ith  radial  position 

Real 

Units 


Meters 


Shooting  Parameters  —  Shooting  parameters  and  related 

variables  computed  by  the  program 


ROOTS (K) 

RATIOA(K) 

RMIN(K) 


kz  for  kth  mode 


a  for  kth  mode 


Final  value  of  minimum  in 
E04JBF  for  kth  mode 


Complex 

Real 

Real 


Meter 


ORTHOG(L, J,K)  Computed  value  of  orthogonal¬ 
ity  between  modes  J  and  K  at 
frequency  L 


Integration  Variables 


The  integration  variables  below 
are  defined  in  Section  3.6 


YVAR(l) 
YVAR( 2 ) 
YVAR( 3 ) 
YVAR(4) 
YPRIME(l) 
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Variable 

Name 


Function 


Units 


IZP_e 


YPRIME(2) 

Y2 

Real 

— 

YPRIME(3) 

Y3 

Real 

— 

YPRIME(4) 

Y4 

Real 

— 

Field  Components  —  Field  values  as 

a  function  of 

r  and 

related  arrays 

used 

in  processing  the  field 

components 

ER(I,L,K) 

Er  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Volts  meter-* 

ET(I,L,K) 

E.  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Volts  meter  * 

EZ( I ,L ,K) 

Ez  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Volts  meter-* 

DR(I,L,K) 

Dr  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Coulombs  meter" 

DT(I,L,K) 

D  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Coulombs  meter" 

DZ(I,L,K) 

D,  at  the 
z 

frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Coulombs  meter" 

BR(I,L,K) 

Br  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Weber  meter 

bt(i,l,k) 

B  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

_2 

Weber  meter 

bz(i,l,k) 

Bz  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Weber  meter 

HR(I,L,K) 

Hr  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Amperes  meter- 

ht(i,l,k) 

H  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Amperes  meter- 

hz(i,l,k) 

Hz  at  the 
frequency 

ith  radial 
L  and  mode 

point 

K 

for 

Real 

Amperes  meter- 
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Variable 

Name 


Function 


Units 


PC(I,L,K) 


PS(I,L,K) 


Value  of  the 
component  Er 
radial  point 
and  mode  K 


Poynting  vector 
*  H  at  the  ith 
for”f requency  L 


Value  of  the  Poynting  vector 
component  E  <  Hr  at  the  ith 
radial  point  for  frequency  L 
and  mode  K 


Tjrpe 


Real  Watts 


Real  Watts 


FLD(I,L,K,N) 

Used  in  processing  the  field 
components i  equivalenced  to  the 
Nth  field  component 

Real 

DRP( I »L,K) 

Derivative  of  Dr  with  respect 
to  r  at  the  ith  radial  point 
for  frequency  L  and  mode  K 

Real 

Coulombs  meter 

ETP(I,L,K) 

Derivative  of  E,  with  respect 
to  r  at  the  ith” radial  point 
for  frequency  L  and  mode  K 

Real 

Volts  meter 

DTP( I ,L,K) 

Derivative  of  D  with  respect 
to  r  at  the  ith* radial  point 
for  frequency  L  and  mode  K 

Real 

Coulombs  meter 

ERDIFF(I,L,K) 

Value  of  the  "jump"  in  Er  at 
the  ith  discontinuity  in  the 
dielectric 

Real 

Volts  meter  1 

DTD1FF( I ,L,K) 

Value  of  the  "jump"  in  D  at 
the  ith  discontinuity  in”the 
dielectric 

Real 

Coulombs  meter 

dzdiff(i,l,k) 

Value  of  the  "jump"  in  Dz  at 
the  ith  discontinuity  in  the 
dielectric 

Real 

Coulombs  meter 

erd(i,l,k) 

"Smoothed"  field  component  Er; 
ERD(I,L,K)  -  ER( I ,L,K)  - 
ERDIFF(I,L,K) 

Real 

Volts  meter  1 

DTD(I,L,K) 

"Smoothed"  field  component  D  ; 
DTD(I,L,K)  -  DT(I,L,K)  -  * 
DTDIFF(I,L,K) 

Real 

Coulombs  meter 

dzd(i,l,k) 

"Smoothed"  field  component  Dz; 
DZD( I ,L,K)  -  DZ( I ,L,K)  - 
DZDIFF(I,L,K) 

Real 

Coulombs  meter" 
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Variable 

Name 


Function 


TyPe 


Units 


D.2  Subroutine  DERIV 


FDRV 

First  derivative  of  the  dielec¬ 

Real 

DREL 

tric  profile  at  a  given  point 

Dielectric  value  at  a  given  point 

Real 

FCTRl 

Normalized 

value  of  r  1 

Real 

— 

FCTR2 

Normalized 

value  of 

z 

Real 

— 

FCTR3 

Normalized 

2  -2 

value  of  a)  c 

Real 

— 

FCTR4 

m^  +  1 

Real 

— 

FCTR5 

DREL-1 

Real 

— 

D.3  Subroutine  PARDRV 

Processing  and  Control  — 

The  variables  FDRV,  DREL,  FCTRl 

,  FCTR2 ,  FCTR3 

PD(1,1) 

PD1,1 

FCTR4,  and  FCTR5  are  used  here  as  in  sub¬ 
routine  DERIV  and  will  not  be  repeated.  The 
following  variables  are  defined  in  Section  3. 

Real  - 

PD(1 ,2) 

PD1 ,2 

Real 

— 

PD( 1 ,3) 

PD1 ,3 

Real 

— 

PD( 1,4) 

PD1 ,4 

Real 

— 

PD( 2,1) 

PD2,1 

Real 

— 

PD(2 ,2) 

PD2  ,2 

Real 

— 

PD( 2,3) 

PD2 ,3 

Real 

— 

PD(2.4) 

PD2,4 

Real 

— 

PD(3 ,1) 

PD3,1 

Real 

— 

PD(3 ,2) 

PD3 ,2 

Real 

— 

PD(3 ,3) 

PD3 ,3 

Real 

— 

PD( 3,4) 

PD3,4 

Real 

— 

tjtf 

*  V 


n 


I 


a 


P 

I 


Var<  able 


Name 

Function 

1X21 

Units 

PD(4 ,1) 

PD4,1 

Real 

— 

PD(4  ,2) 

PD4,2 

Real 

— 

PD(4 ,3) 

PD4,3 

Real 

— 

PD(4  ,4) 

PD4,4 

Real 

— 

D.4  Function 

FNCT1 

XEND 

Normalized  endpoint 

used  by  DGEAR 

Real 

— 

FNCT1 

Functional  value 
of  the  following 

equal  to  one 
equations : 

Real 

— 

™0  ,n 

and  TEM  modes 

:  equation 

93 

TE0,n 

modes 

:  equation 

94 

-  ■ 

D.5  Function  FNCT2 

XEND  Normalized  endpoint  used  by  DGEAR  Real  - 

FC  Functional  value  equal  to  Eq.  95  Peal  - 

D.6  Subroutine  NORMAL 

Spline  Arrays  —  The  following  arrays  store  the  computed  spline  coeffi¬ 
cients  generated  by  the  routine  ICSCCU  and  evaluated  by 
ICSEVU  and  DCSEVU.  The  arrays  below  are  also  used  in 
subroutine  ORTHO  and  function  CX1NT. 


ERS( 1,3,2) 

Spline  coefficients 
ith  radial  position 

for 

Er 

at 

Real 

ETS( I ,3,2) 

Spline  coefficients 
ith  radial  position 

for 

at 

Real 

— 

HRS( 1,3,2) 

Spline  coefficients 
ith  radial  position 

for 

Hr 

at 

Real 

— 

HTS( 1,3,2) 


Spline  coefficients  for  H  at 
ith  radial  position 


Real 


0'VV’*Vv. 


Variable 

Name 


Function 


D.7  IMSL  and  NAG  Routines 


DGEAR 


Number  of  first  order  differ¬ 
ential  equations 


Integer 


ZREAL1 


On  input,  XP  supplies  the  ini-  Real 
tial  value  and  is  used  only  on 
the  first  call.  On  output,  XP 
is  replaced  with  the  current 
value  of  the  independent  vari¬ 
able  at  which  integration  has 
been  completed 

Integer  work  array  of  length  N  Integer 

Work  array  of  length  4*N+METHT+  Real 
MITERT 


Spread  criteria  for  multiple 
roots 


E04JBE 


Used  to  restart  a  computation 
when  multiple  roots  are  desired 


Number  of  variables 


INPRINT  Specifies  the  frequency  with 

which  subroutine  MONIT  is  to 
be  called.  There  are  three 
options : 


Integer 


IPRINT  >  0:  MONIT  called  once 
every  IPRINT  iter¬ 
ations 


IPRINT  -  0: 


MONIT  called  at 
final  point  only 


IPRINT  <  0: 


MONIT  not  called 
at  all 


Variable 

Name 


Function 


Units 


II 21  _ 


LOCSH 


Specifies  whether  or  not  the 
user  wishes  a  "local  search" 
(TRUE  or  FALSE,  respectively) 
be  performed  when  a  point  is 
found  which  is  thought  to  be 
a  constrained  minimum 


Logical 


Upon  judicious  choice  for  the 
range  0.0  <  ETA  <  1.0,  the 
linear  minimization  process 
is  made  more  efficient 


m 


I BOUND 


Specifies  whether  the  problem 
is  unconstrained  or  bounded. 
The  options  are: 


Integer 


I BOUND  *  0:  Variables  are 

bounded  and  the 
user  supplies  L. 
and  Uj  ^ 

I BOUND  *  1:  Problem  is  uncon¬ 
strained  and  the 
function  FUNCT  is 
called  N  times 


I BOUND  -  2:  Variables  are 
bounded 


Specifies  the  actual  dimension 
of  the  NAG  subroutine  HESL 
called  by  E04JBF 


Integer 


m 

m 


Workspace  array  Integer 

Specifies  the  actual  dimension  Integer 
of  IW  as  declared  in  the  (sub) 
program  from  which  E04JBF  is 
called 


33 

wd 

Sfi 

«'{»$• 

m 

‘Va 


Workspace  array  of  at  least  9*N 

Specifies  the  actual  dimension 
of  W  as  declared  in  the  (sub) 
program  from  which  E04JBF  is 
called 


Integer 
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Si  V, 


Variable 

Name  Function 


DCADRE 


Desired  relative  error  in  the 
answer 


RERR 


u\ 

» 


The 

values  of 

APPENDIX  E 

VALUES  OF  THE  RADIAL 

the  radial  po*-  ts  at 

POINTS 

which  the 

dielectric  profile 

was  defined 

for  both 

cases  in  Chapter  5  are  presented  below  (in  centi- 

meters) : 

rl 

m 

1.0000 

r24 

m 

1.5000 

r2 

m 

1.0400 

r25 

m 

1.5010 

r3 

m 

1.0800 

r26 

m 

1.5050 

r4 

m 

1.1200 

r27 

m 

1.5075 

r5 

m 

1.1600 

r28 

m 

1.5100 

r6 

m 

1.2000 

r29 

m 

1.5150 

r7 

m 

1.2400 

r30 

m 

1.5175 

r8 

m 

1.2800 

r31 

m 

1.5200 

r9 

m 

1.3200 

r32 

m 

1.5250 

r10 

s 

1.3600 

r33 

m 

1.5300 

rll 

m 

1.4000 

r34 

m 

1.5400 

r12 

m 

1.4400 

r35 

m 

1.5600 

r13 

m 

1.4600 

r36 

m 

1.6000 

r14 

m 

1.4650 

r37 

m 

1.6400 

r15 

m 

1.4700 

r38 

m 

1.6800 

r16 

m 

1.4750 

r39 

m 

1.7200 

r17 

m 

1.4800 

r40 

m 

1.7600 

r18 

m 

1.4850 

r41 

m 

1.8000 

r19 

m 

1.4900 

r42 

m 

1.8400 

o 

CM 

U 

m 

1.4925 

r43 

m 

1.8800 

r21 

m 

1.4950 

r44 

m 

1.9200 

r22 

m 

1.4990 

r45 

m 

1.9600 

r-ii 

m 

1.5000 

r/,t 

m 

2.0000 

% 
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